HarrisExample.cxx 6.31 KB
 Sébastien Dinot committed Mar 08, 2017 1 /*  Julien Michel committed Jan 14, 2019 2  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)  Sébastien Dinot committed Mar 08, 2017 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 committed Apr 04, 2006 20   Patrick Imbo committed Jun 29, 2006 21   Guillaume Pasero committed Mar 20, 2019 22 #include "otbMacro.h"  Jordi Inglada committed May 18, 2006 23 #include "otbImage.h"  Jordi Inglada committed Apr 04, 2006 24 25 26 27  #include "otbImageFileReader.h" #include "otbImageFileWriter.h"  Victor Poughon committed Feb 20, 2019 28 29 30 31 /* Example usage: ./HarrisExample Input/ROISpot5.png Output/ROISpot5Harris.png 1.5 2 0.1 */  Jordi Inglada committed Apr 04, 2006 32   Jordi Inglada committed Jul 03, 2006 33 // This example illustrates the use of the \doxygen{otb}{HarrisImageFilter}.  Jordi Inglada committed Apr 04, 2006 34 //  Emmanuel Christophe committed Dec 06, 2008 35 // The first step required to use this filter is to include its header file.  Jordi Inglada committed Apr 04, 2006 36   Jordi Inglada committed Apr 13, 2006 37 #include "otbHarrisImageToPointSetFilter.h"  Jordi Inglada committed Apr 04, 2006 38 39 #include "itkRescaleIntensityImageFilter.h"  Victor Poughon committed Feb 20, 2019 40 int main(int argc, char* argv[])  Jordi Inglada committed Apr 04, 2006 41 {  OTB Bot committed Apr 01, 2010 42  if (argc != 6)  Victor Poughon committed Feb 20, 2019 43  {  Jordi Inglada committed Apr 04, 2006 44  std::cerr << "Usage: " << argv[0] << " inputImageFile ";  Emmanuel Christophe committed Dec 06, 2008 45  std::cerr << " outputHarrisImageFile sigmaD sigmaI alpha" << std::endl;  Jordi Inglada committed Apr 04, 2006 46  return EXIT_FAILURE;  Victor Poughon committed Feb 20, 2019 47  }  Emmanuel Christophe committed Jan 31, 2009 48   Victor Poughon committed Feb 20, 2019 49 50  const char* inputFilename = argv[1]; const char* outputFilename = argv[2];  Emmanuel Christophe committed Jan 31, 2009 51   Victor Poughon committed Feb 20, 2019 52 53 54  double SigmaD((double)::atof(argv[3])); double SigmaI((double)::atof(argv[4])); double Alpha((double)::atof(argv[5]));  Emmanuel Christophe committed Jan 31, 2009 55   Victor Poughon committed May 16, 2019 56 57 58  using InputPixelType = float; const unsigned int Dimension = 2; using OutputPixelType = unsigned char;  Emmanuel Christophe committed Jan 31, 2009 59   Victor Poughon committed May 16, 2019 60 61  using InputImageType = otb::Image; using OutputImageType = otb::Image;  Emmanuel Christophe committed Jan 31, 2009 62   Victor Poughon committed May 16, 2019 63  using ReaderType = otb::ImageFileReader;  Emmanuel Christophe committed Jan 31, 2009 64 65 66 67 68  // The \doxygen{otb}{HarrisImageFilter} is templated over the // input and output image types, so we start by // defining:  Victor Poughon committed May 16, 2019 69 70  using HarrisFilterType = otb::HarrisImageFilter; using RescalerType = itk::RescaleIntensityImageFilter;  Emmanuel Christophe committed Jan 31, 2009 71   Victor Poughon committed May 16, 2019 72  using WriterType = otb::ImageFileWriter;  Emmanuel Christophe committed Jan 31, 2009 73   OTB Bot committed Apr 01, 2010 74 75  ReaderType::Pointer reader = ReaderType::New(); WriterType::Pointer writer = WriterType::New();  Victor Poughon committed Feb 20, 2019 76  HarrisFilterType::Pointer harris = HarrisFilterType::New();  OTB Bot committed Apr 01, 2010 77  RescalerType::Pointer rescaler = RescalerType::New();  Emmanuel Christophe committed Jan 31, 2009 78   OTB Bot committed Apr 01, 2010 79 80  reader->SetFileName(inputFilename); writer->SetFileName(outputFilename);  Emmanuel Christophe committed Jan 31, 2009 81   OTB Bot committed Apr 01, 2010 82  harris->SetInput(reader->GetOutput());  Emmanuel Christophe committed Jan 31, 2009 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: // // \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)&  Julien Malik committed Jun 28, 2012 95 96  // L_y^2(\mathbf{x},\sigma_D) \end{array}\right] //  Emmanuel Christophe committed Jan 31, 2009 97 98  // The output of the detector is $$det(\mu) - \alpha trace^2(\mu).$$  OTB Bot committed Apr 01, 2010 99 100 101  harris->SetSigmaD(SigmaD); harris->SetSigmaI(SigmaI); harris->SetAlpha(Alpha);  Emmanuel Christophe committed Jan 31, 2009 102   OTB Bot committed Apr 01, 2010 103 104  rescaler->SetOutputMinimum(itk::NumericTraits::min()); rescaler->SetOutputMaximum(itk::NumericTraits::max());  Emmanuel Christophe committed Jan 31, 2009 105   OTB Bot committed Apr 01, 2010 106 107  rescaler->SetInput(harris->GetOutput()); writer->SetInput(rescaler->GetOutput());  Emmanuel Christophe committed Jan 31, 2009 108 109  writer->Update();  Jordi Inglada committed Apr 04, 2006 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 committed May 14, 2008 116  // \itkcaption[Harris Filter Application]{Result of applying the  Emmanuel Christophe committed Dec 06, 2008 117  // \doxygen{otb}{HarrisImageFilter} to a Spot 5 image.}  Jordi Inglada committed Apr 04, 2006 118 119 120  // \label{fig:Harris} // \end{figure} //  Jordi Inglada committed Jul 03, 2006 121  // The output of the \doxygen{otb}{HarrisImageFilter} is an image  OTB Bot committed Jan 31, 2011 122  // where, for each pixel, we obtain the intensity of the  Jordi Inglada committed Apr 13, 2006 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  Jordi Inglada committed Jul 03, 2006 126  // \doxygen{otb}{HarrisImageToPointSetFilter}. This filter is only  Jordi Inglada committed Apr 13, 2006 127  // templated over the input image type, the output being a  Jordi Inglada committed Jul 03, 2006 128  // \doxygen{itk}{PointSet} with pixel type equal to the image pixel type.  Jordi Inglada committed Apr 04, 2006 129   Victor Poughon committed May 16, 2019 130  using FunctionType = otb::HarrisImageToPointSetFilter;  Emmanuel Christophe committed Jan 31, 2009 131 132  // We declare now the filter and a pointer to the output point set.  Victor Poughon committed May 16, 2019 133  using OutputPointSetType = FunctionType::OutputPointSetType;  Emmanuel Christophe committed Jan 31, 2009 134   Victor Poughon committed Feb 20, 2019 135 136  FunctionType::Pointer harrisPoints = FunctionType::New(); OutputPointSetType::Pointer pointSet = OutputPointSetType::New();  Emmanuel Christophe committed Jan 31, 2009 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 committed Apr 01, 2010 142 143 144 145 146  harrisPoints->SetInput(0, reader->GetOutput()); harrisPoints->SetSigmaD(SigmaD); harrisPoints->SetSigmaI(SigmaI); harrisPoints->SetAlpha(Alpha); harrisPoints->SetLowerThreshold(10);  Emmanuel Christophe committed Jan 31, 2009 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.  Victor Poughon committed May 16, 2019 158 159 160 161  using ContainerType = OutputPointSetType::PointsContainer; ContainerType* pointsContainer = pointSet->GetPoints(); using IteratorType = ContainerType::Iterator; IteratorType itList = pointsContainer->Begin();  Emmanuel Christophe committed Jan 31, 2009 162 163 164  // And we get the points coordinates  OTB Bot committed Apr 01, 2010 165  while (itList != pointsContainer->End())  Victor Poughon committed Feb 20, 2019 166  {  Victor Poughon committed May 16, 2019 167 168  using OutputPointType = OutputPointSetType::PointType; OutputPointType pCoordinate = (itList.Value());  Guillaume Pasero committed Mar 20, 2019 169  otbLogMacro(Debug, << pCoordinate);  Emmanuel Christophe committed Jan 31, 2009 170  ++itList;  Victor Poughon committed Feb 20, 2019 171  }  Jordi Inglada committed Apr 13, 2006 172   Jordi Inglada committed Apr 04, 2006 173 174  return EXIT_SUCCESS; }