AVIMultiChannelRAndGAndNIRVegetationIndexImageFilter.cxx 8.19 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.
 */
20

21

22 23 24 25 26 27 28
/* Example usage:
./AVIMultiChannelRAndGAndNIRVegetationIndexImageFilter Input/verySmallFSATSW.tif \
                                                       Output/AVIMultiChannelRAndGAndNIRVegetationIndex.tif \
                                                       Output/pretty_FSATSW.png \
                                                       Output/pretty_AVIMultiChannelRAndGAndNIRVegetationIndex.png \
                                                       3 \
                                                       2 \
Julien Michel's avatar
Julien Michel committed
29
                                                       4
30
*/
31 32 33 34 35 36 37


// \index{otb::VegetationIndex}
// \index{otb::VegetationIndex!header}
//
//
// The following example illustrates the use of the
38
// itk::UnaryFunctorImageFilter with the
39
// use of the Angular Vegetation Index (AVI).
40
// The equation for the Angular Vegetation Index involves the gren, red
Jordi Inglada's avatar
Jordi Inglada committed
41
// and near infra-red bands. $\lambda_1$, $\lambda_2$ and $\lambda_3$ are the mid-band
42
// wavelengths for the green, red and NIR bands and $\tan^{-1}$ is the arctangent function.
43 44 45 46
//
// The AVI expression is
//
// \begin{equation}
Jordi Inglada's avatar
Jordi Inglada committed
47
// \mathbf{A_1} = \frac{\lambda_3-\lambda_2}{\lambda_2}
48 49
// \end{equation}
// \begin{equation}
Jordi Inglada's avatar
Jordi Inglada committed
50
// \mathbf{A_2} = \frac{\lambda_2-\lambda_1}{\lambda_2}
51 52 53
// \end{equation}
//
// \begin{equation}
Jordi Inglada's avatar
Jordi Inglada committed
54
// \mathbf{AVI} = \tan^{-1}\left({\frac{A_1}{NIR-R}}\right) + \tan^{-1}\left({\frac{A_2}{G-R}}\right)
55 56 57 58 59
// \end{equation}
//
// For more details, refer to Plummer work \cite{AVI}.
//
//
60
// Let's look at the minimal code required to use this
61
// algorithm.
62

63 64
#include "otbVegetationIndicesFunctor.h"
#include "itkUnaryFunctorImageFilter.h"
65 66 67 68 69 70 71 72 73

#include "otbImage.h"
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "otbVectorRescaleIntensityImageFilter.h"
#include "otbMultiChannelExtractROI.h"
#include "itkThresholdImageFilter.h"

74
int main(int argc, char* argv[])
75
{
Julien Michel's avatar
Julien Michel committed
76
  if (argc < 8)
77
  {
78 79
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
Julien Michel's avatar
Julien Michel committed
80
    std::cerr << " inputImage , outputImage , prettyInput , prettyOutput , redChannel , greenChannel , nirChannel ," << std::endl;
81
    return 1;
82
  }
83 84 85 86 87

  // The image types are now defined using pixel types and
  // dimension. The input image is defined as an \doxygen{otb}{VectorImage},
  // the output is a \doxygen{otb}{Image}.

88 89 90 91 92
  const unsigned int Dimension = 2;
  using InputPixelType         = double;
  using OutputPixelType        = float;
  using InputImageType         = otb::VectorImage<InputPixelType, Dimension>;
  using OutputImageType        = otb::Image<OutputPixelType, Dimension>;
93 94

  // We instantiate reader and writer types
95 96
  using ReaderType = otb::ImageFileReader<InputImageType>;
  using WriterType = otb::ImageFileWriter<OutputImageType>;
97 98 99 100

  // The AVI (Angular Vegetation Index) is
  // instantiated using the image pixel types as template parameters.

101
  using FunctorType = otb::Functor::AVI<InputPixelType, OutputPixelType>;
102 103

  // The
104
  // \doxygen{itk}{UnaryFunctorImageFilter}
105 106 107
  // type is defined using the image types and the AVI functor as
  // template parameters. We then instantiate the filter itself.

108
  using AVIImageFilterTypeType = itk::UnaryFunctorImageFilter<InputImageType, OutputImageType, FunctorType>;
109

110
  AVIImageFilterTypeType::Pointer filter = AVIImageFilterTypeType::New();
111 112 113 114 115 116

  ReaderType::Pointer reader = ReaderType::New();
  WriterType::Pointer writer = WriterType::New();

  //  Now the input image is set and a name is given to the output image.

OTB Bot's avatar
OTB Bot committed
117 118
  reader->SetFileName(argv[1]);
  writer->SetFileName(argv[2]);
119 120 121

  // The three used index bands (red, green and NIR) are declared.

Julien Michel's avatar
Julien Michel committed
122 123 124
  filter->GetFunctor().SetBandIndex(CommonBandNames::RED, ::atoi(argv[5]));
  filter->GetFunctor().SetBandIndex(CommonBandNames::GREEN, ::atoi(argv[6]));
  filter->GetFunctor().SetBandIndex(CommonBandNames::NIR, ::atoi(argv[7]));
125 126 127 128

  // The filter input is linked to the reader output and
  // the filter output is linked to the writer input.

OTB Bot's avatar
OTB Bot committed
129
  filter->SetInput(reader->GetOutput());
130

OTB Bot's avatar
OTB Bot committed
131
  writer->SetInput(filter->GetOutput());
132 133 134 135 136 137

  //  The invocation of the \code{Update()} method on the writer triggers the
  //  execution of the pipeline.  It is recommended to place update calls in a
  //  \code{try/catch} block in case errors occur and exceptions are thrown.

  try
138
  {
139
    writer->Update();
140
  }
OTB Bot's avatar
OTB Bot committed
141
  catch (itk::ExceptionObject& excep)
142
  {
143 144
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
145
  }
OTB Bot's avatar
OTB Bot committed
146
  catch (...)
147
  {
148 149
    std::cout << "Unknown exception !" << std::endl;
    return EXIT_FAILURE;
150
  }
151 152

  // Pretty image creation for the printing
153 154 155 156 157 158 159
  using OutputPrettyImageType       = otb::Image<unsigned char, Dimension>;
  using OutputVectorPrettyImageType = otb::VectorImage<unsigned char, Dimension>;
  using WriterVectorPrettyType      = otb::ImageFileWriter<OutputVectorPrettyImageType>;
  using WriterPrettyType            = otb::ImageFileWriter<OutputPrettyImageType>;
  using RescalerType                = itk::RescaleIntensityImageFilter<OutputImageType, OutputPrettyImageType>;
  using VectorRescalerType          = otb::VectorRescaleIntensityImageFilter<InputImageType, OutputVectorPrettyImageType>;
  using ChannelExtractorType        = otb::MultiChannelExtractROI<unsigned char, unsigned char>;
160 161 162 163

  VectorRescalerType::Pointer     vectRescaler     = VectorRescalerType::New();
  ChannelExtractorType::Pointer   selecter         = ChannelExtractorType::New();
  WriterVectorPrettyType::Pointer vectPrettyWriter = WriterVectorPrettyType::New();
OTB Bot's avatar
OTB Bot committed
164 165

  OutputVectorPrettyImageType::PixelType minimum, maximum;
166 167 168 169 170 171
  minimum.SetSize(reader->GetOutput()->GetNumberOfComponentsPerPixel());
  maximum.SetSize(reader->GetOutput()->GetNumberOfComponentsPerPixel());
  minimum.Fill(0);
  maximum.Fill(255);
  vectRescaler->SetOutputMinimum(minimum);
  vectRescaler->SetOutputMaximum(maximum);
172
  //  vectRescaler->SetClampThreshold(1);
OTB Bot's avatar
OTB Bot committed
173
  vectRescaler->SetInput(reader->GetOutput());
174 175 176

  selecter->SetInput(vectRescaler->GetOutput());
  selecter->SetChannel(3);
Jordi Inglada's avatar
Jordi Inglada committed
177 178
  selecter->SetChannel(2);
  selecter->SetChannel(1);
179

OTB Bot's avatar
OTB Bot committed
180 181
  vectPrettyWriter->SetFileName(argv[3]);
  vectPrettyWriter->SetInput(selecter->GetOutput());
182

183
  using ThresholderType = itk::ThresholdImageFilter<OutputImageType>;
184 185

  ThresholderType::Pointer thresholder = ThresholderType::New();
OTB Bot's avatar
OTB Bot committed
186 187 188
  thresholder->SetInput(filter->GetOutput());
  thresholder->SetOutsideValue(1.0);
  thresholder->ThresholdOutside(-1.0, 0.05);
189 190 191 192
  thresholder->Update();

  RescalerType::Pointer     rescaler     = RescalerType::New();
  WriterPrettyType::Pointer prettyWriter = WriterPrettyType::New();
OTB Bot's avatar
OTB Bot committed
193
  rescaler->SetInput(thresholder->GetOutput());
194 195
  rescaler->SetOutputMinimum(0);
  rescaler->SetOutputMaximum(255);
OTB Bot's avatar
OTB Bot committed
196 197
  prettyWriter->SetFileName(argv[4]);
  prettyWriter->SetInput(rescaler->GetOutput());
198 199

  try
200
  {
201 202
    prettyWriter->Update();
    vectPrettyWriter->Update();
203
  }
OTB Bot's avatar
OTB Bot committed
204
  catch (itk::ExceptionObject& excep)
205
  {
206 207
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
208
  }
OTB Bot's avatar
OTB Bot committed
209
  catch (...)
210
  {
211 212
    std::cout << "Unknown exception !" << std::endl;
    return EXIT_FAILURE;
213
  }
214 215 216 217

  return EXIT_SUCCESS;

  // Let's now run this example using as input the image
218
  // \code{verySmallFSATSW.tif} provided in the
219 220 221 222
  // directory \code{Examples/Data}.
  //
  //
  // \begin{figure} \center
223
  // \includegraphics[width=0.24\textwidth]{pretty_FSATSW.eps}
224 225
  // \includegraphics[width=0.24\textwidth]{pretty_AVIMultiChannelRAndGAndNIRVegetationIndex.eps}
  // \itkcaption[AVI Example]{AVI result on the right with the left image in input.}
226
  // \label{fig:AVIMultiChannelRAndGAndNIRIndexImageFilter}
227 228
  // \end{figure}
}