HarrisExample.cxx 7.95 KB
Newer Older
Jordi Inglada's avatar
Jordi Inglada committed
1 2
/*=========================================================================

Patrick Imbo's avatar
nomsg  
Patrick Imbo committed
3 4 5 6 7 8 9 10 11 12 13 14 15
  Program:   ORFEO Toolbox
  Language:  C++
  Date:      $Date$
  Version:   $Revision$


  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
  See OTBCopyright.txt for details.


     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.
Jordi Inglada's avatar
Jordi Inglada committed
16 17

=========================================================================*/
Patrick Imbo's avatar
nomsg  
Patrick Imbo committed
18

Jordi Inglada's avatar
Jordi Inglada committed
19 20 21 22 23
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include "itkExceptionObject.h"
24
#include "otbImage.h"
Jordi Inglada's avatar
Jordi Inglada committed
25 26 27 28 29 30 31 32 33 34 35 36

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

//  Software Guide : BeginCommandLineArgs
//    INPUTS: {ROISpot5.png}
//    OUTPUTS: {ROISpot5Harris.png}
//    1.5 2 0.1
//  Software Guide : EndCommandLineArgs

// Software Guide : BeginLatex
//
37
// This example illustrates the use of the \doxygen{otb}{HarrisImageFilter}.
Jordi Inglada's avatar
Jordi Inglada committed
38 39 40 41 42 43 44 45
// 
// The first step required to use this filter is to include its header file. 
//
// Software Guide : EndLatex 

// Software Guide : BeginCodeSnippet
#include "otbHarrisImageFilter.h"
// Software Guide : EndCodeSnippet
46
#include "otbHarrisImageToPointSetFilter.h"
Jordi Inglada's avatar
Jordi Inglada committed
47 48
#include "itkRescaleIntensityImageFilter.h"

49
int main(int argc, char *argv[] )
Jordi Inglada's avatar
Jordi Inglada committed
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
{
    if( argc != 6 )
    {
    std::cerr << "Usage: " << argv[0] << " inputImageFile ";
    std::cerr << " outputHarrisImageFile sigmaD sigmaI alpha" << std::endl;  
    return EXIT_FAILURE;
    }

    const char * inputFilename  = argv[1];
    const char * outputFilename = argv[2];
        
    double SigmaD((double)::atof(argv[3]));
    double SigmaI((double)::atof(argv[4]));
    double Alpha((double)::atof(argv[5]));
	        
    typedef float                                   InputPixelType;
    const   unsigned int        	                        Dimension = 2;
    typedef unsigned char   	                        OutputPixelType;
	
69 70
    typedef otb::Image< InputPixelType,  Dimension >   InputImageType;
    typedef otb::Image< OutputPixelType, Dimension >   OutputImageType;
Jordi Inglada's avatar
Jordi Inglada committed
71 72 73 74 75

    typedef otb::ImageFileReader< InputImageType  >    ReaderType;

    //  Software Guide : BeginLatex
    //
76
    //  The \doxygen{otb}{HarrisImageFilter} is templated over the
Jordi Inglada's avatar
Jordi Inglada committed
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104
    //  input and output image types, so we start by
    //  defining:
    //
    //  Software Guide : EndLatex 
    
    // Software Guide : BeginCodeSnippet
    typedef otb::HarrisImageFilter<InputImageType,
                                     InputImageType>   HarrisFilterType;

    // Software Guide : EndCodeSnippet
    typedef itk::RescaleIntensityImageFilter
                   < InputImageType, OutputImageType > RescalerType; 

    typedef otb::ImageFileWriter< OutputImageType >    WriterType;
  

    ReaderType::Pointer reader   = ReaderType::New();
    WriterType::Pointer writer   = WriterType::New();
    HarrisFilterType::Pointer harris = HarrisFilterType::New();
    RescalerType::Pointer rescaler = RescalerType::New();
	
    reader->SetFileName( inputFilename  );
    writer->SetFileName( outputFilename );

    harris->SetInput( reader->GetOutput() );
    
    //  Software Guide : BeginLatex
    //
105
    // The \doxygen{otb}{HarrisImageFilter} needs some parameters to
Jordi Inglada's avatar
Jordi Inglada committed
106 107 108 109 110
    // 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
Patrick Imbo's avatar
Patrick Imbo committed
111
    // scale). This allows the computation of the following matrix:
Jordi Inglada's avatar
Jordi Inglada committed
112 113 114
    // \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) &
Jordi Inglada's avatar
Jordi Inglada committed
115
    // L_xL_y(\mathbf{x},\sigma_D)\\ L_xL_y(\mathbf{x},\sigma_D)&
Jordi Inglada's avatar
Jordi Inglada committed
116
    // L_y^2(\mathbf{x},\sigma_D) \end{array}\right] \end{equation}
Jordi Inglada's avatar
Jordi Inglada committed
117
    // The output of the detector is $$det(\mu) - \alpha trace^2(\mu).$$
Jordi Inglada's avatar
Jordi Inglada committed
118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
    //
    //  Software Guide : EndLatex

    // Software Guide : BeginCodeSnippet
    harris->SetSigmaD( SigmaD );
    harris->SetSigmaI( SigmaI );
    harris->SetAlpha( Alpha );
    // Software Guide : EndCodeSnippet
    
    rescaler->SetOutputMinimum( itk::NumericTraits< OutputPixelType >::min());
    rescaler->SetOutputMaximum( itk::NumericTraits< OutputPixelType >::max());

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

    //  Software Guide : BeginLatex
  // 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
141
  // \itkcaption[Harris Filter Application]{Result of applying the
142
  // \doxygen{otb}{HarrisImageFilter} to a Spot 5 image.} 
Jordi Inglada's avatar
Jordi Inglada committed
143 144 145
  // \label{fig:Harris}
  // \end{figure}
  //
146
  // The output of the \doxygen{otb}{HarrisImageFilter} is an image
147 148 149 150
  // where,for each pixel, we obtain the intensity of the
  // 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
151
  // \doxygen{otb}{HarrisImageToPointSetFilter}. This filter is only
152
  // templated over the input image type, the output being a
153
  // \doxygen{itk}{PointSet} with pixel type equal to the image pixel type.
154
  // 
Jordi Inglada's avatar
Jordi Inglada committed
155 156
  //  Software Guide : EndLatex 

157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174
    // Software Guide : BeginCodeSnippet
    typedef otb::HarrisImageToPointSetFilter<InputImageType>   FunctionType;
    // Software Guide : EndCodeSnippet

    //  Software Guide : BeginLatex
    //
    //  We declare now the filter and a pointer to the output point set.
    //  Software Guide : EndLatex 
    // Software Guide : BeginCodeSnippet
    typedef FunctionType::OutputPointSetType   OutputPointSetType;


    FunctionType::Pointer        harrisPoints    = FunctionType::New();
    OutputPointSetType::Pointer  pointSet = OutputPointSetType::New();
    // Software Guide : EndCodeSnippet

    //  Software Guide : BeginLatex
    //
175 176
    //  The \doxygen{otb}{HarrisImageToPointSetFilter} takes the same
    // parameters as the \doxygen{otb}{HarrisImageFilter} and an
177 178 179 180 181 182 183 184 185 186
    // additional parameter : the threshold for the point selection.
    //
    //  Software Guide : EndLatex
    
    // Software Guide : BeginCodeSnippet

    harrisPoints->SetInput( 0,reader->GetOutput() );
    harrisPoints->SetSigmaD( SigmaD );
    harrisPoints->SetSigmaI( SigmaI );
    harrisPoints->SetAlpha( Alpha );
187
    harrisPoints->SetLowerThreshold( 10 );
188 189 190 191 192 193 194 195 196 197 198 199
    pointSet = harrisPoints->GetOutput();
    // Software Guide : EndCodeSnippet

    harrisPoints->Update();


    //  Software Guide : BeginLatex
    //
    //  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
200
    //  information on using \doxygen{itk}{PointSet}s) and declaring
201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233
    //  an iterator to it.
    //
    //  Software Guide : EndLatex
    
    // Software Guide : BeginCodeSnippet

    typedef OutputPointSetType::PointsContainer ContainerType;
    ContainerType* pointsContainer = pointSet->GetPoints();
    typedef ContainerType::Iterator IteratorType;
    IteratorType itList = pointsContainer->Begin();
    // Software Guide : EndCodeSnippet

    //  Software Guide : BeginLatex
    //
    //  And we get the points coordinates 
    //
    //  Software Guide : EndLatex
    
    // Software Guide : BeginCodeSnippet


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

    // Software Guide : EndCodeSnippet



Jordi Inglada's avatar
Jordi Inglada committed
234 235 236 237

  return EXIT_SUCCESS;
}