ARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter.cxx 9.56 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 29 30 31
/* Example usage:
./ARVIMultiChannelRAndBAndNIRVegetationIndexImageFilter Input/VegetationIndex.hd \
                                                        Output/ARVIMultiChannelRAndBAndNIRVegetationIndex.tif \
                                                        Output/pretty_VegetationIndex.png \
                                                        Output/pretty_ARVIMultiChannelRAndBAndNIRVegetationIndex.png \
                                                        1 \
                                                        3 \
                                                        2 \
                                                        0.6
*/
32 33


34 35
// \index{otb::VegetationIndices}
// \index{otb::VegetationIndices!header}
36 37
//
//
38
// The following example illustrates the use of the
39
//  \doxygen{itk}{UnaryFunctorImageFilter} with the
40
// use of the Atmospherically Resistant Vegetation Index (ARVI) \subdoxygen{otb}{Functor}{ARVI}.  ARVI
Julien Malik's avatar
Julien Malik committed
41
// is an improved version of the NDVI that is more robust to the
42 43 44 45 46 47 48 49 50
// atmospheric effect.  In addition to the red and NIR channels (used
// in the NDVI), the ARVI takes advantage of the presence of the blue
// channel to accomplish a self-correction process for the atmospheric
// effect on the red channel. For this, it uses the difference in the
// radiance between the blue and the red channels to correct the
// radiance in the red channel.  Let's define $\rho_{NIR}^{*}$,
// $\rho_{r}^{*}$, $\rho_{b}^{*}$ the normalized radiances (that is to
// say the radiance normalized to reflectance units) of red, blue and
// NIR channels respectively.  $\rho_{rb}^{*}$ is defined as
51
//
52 53 54 55 56 57 58 59 60 61 62
// \begin{equation}
// \rho_{rb}^{*} = \rho_{r}^{*} - \gamma*(\rho_{b}^{*} - \rho_{r}^{*})
// \end{equation}
//
// The ARVI expression is
//
// \begin{equation}
// \mathbf{ARVI} = \frac{\rho_{NIR}^{*}-\rho_{rb}^{*}}{\rho_{NIR}^{*}+\rho_{rb}^{*}}
// \end{equation}
//
//
63
// This formula can be simplified with :
64 65 66 67 68
//
// \begin{equation}
// \mathbf{ARVI} = \frac{ L_{NIR}-L_{rb} }{ L_{NIR}+L_{rb} }
// \end{equation}
//
69
// For more details, refer to Kaufman and Tanre' work \cite{ARVI}.
70
//
71 72 73 74 75 76
//  \relatedClasses
//  \begin{itemize}
//  \item \subdoxygen{otb}{Functor}{TSARVI}
//  \item \subdoxygen{otb}{Functor}{EVI}
//  \end{itemize}

77 78
#include "itkUnaryFunctorImageFilter.h"
#include "otbVegetationIndicesFunctor.h"
79 80 81 82 83 84 85

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

88
int main(int argc, char* argv[])
89
{
OTB Bot's avatar
STYLE  
OTB Bot committed
90
  if (argc < 8)
91
  {
92 93
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
94
    std::cerr << " inputImage , outputImage , prettyInput , prettyOutput , redChannel , blueChannel , nirChannel , gama" << std::endl;
95
    return 1;
96
  }
97 98 99

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

102
  const unsigned int                                  Dimension = 2;
103 104
  typedef double                                      InputPixelType;
  typedef float                                       OutputPixelType;
OTB Bot's avatar
STYLE  
OTB Bot committed
105 106
  typedef otb::VectorImage<InputPixelType, Dimension> InputImageType;
  typedef otb::Image<OutputPixelType, Dimension>      OutputImageType;
107

108 109
  // We instantiate reader and writer types
  typedef otb::ImageFileReader<InputImageType>  ReaderType;
110
  typedef otb::ImageFileWriter<OutputImageType> WriterType;
111

112 113
  // The ARVI (Atmospherically Resistant Vegetation Index) is
  // instantiated using the image pixel types as template parameters.
Emmanuel Christophe's avatar
Emmanuel Christophe committed
114
  // Note that we also can use other functors which operate with the
115
  // Red, Blue and Nir channels such as EVI, ARVI and TSARVI.
116

117
  typedef otb::Functor::ARVI<InputPixelType, InputPixelType, InputPixelType, OutputPixelType> FunctorType;
118

119
  // The
120
  // \doxygen{itk}{UnaryFunctorImageFilter}
121 122
  // type is defined using the image types and the ARVI functor as
  // template parameters. We then instantiate the filter itself.
123

124
  typedef itk::UnaryFunctorImageFilter<InputImageType, OutputImageType, FunctorType> ArviImageFilterType;
125

126
  ArviImageFilterType::Pointer filter = ArviImageFilterType::New();
127

128 129 130 131 132
  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
STYLE  
OTB Bot committed
133 134
  reader->SetFileName(argv[1]);
  writer->SetFileName(argv[2]);
135

136 137
  // The three used index bands (red, blue and NIR) are declared.

138 139 140
  filter->GetFunctor().SetRedIndex(::atoi(argv[5]));
  filter->GetFunctor().SetBlueIndex(::atoi(argv[6]));
  filter->GetFunctor().SetNIRIndex(::atoi(argv[7]));
141

142
  // The $\gamma$ parameter is set. The
143
  // \doxygen{otb::Functor}{ARVI}
144 145
  // class sets the default value of $\gamma$ to $0.5$.  This parameter
  // is used to reduce the atmospheric effect on a global scale.
146 147

  filter->GetFunctor().SetGamma(::atof(argv[8]));
148

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

OTB Bot's avatar
STYLE  
OTB Bot committed
152
  filter->SetInput(reader->GetOutput());
153

OTB Bot's avatar
STYLE  
OTB Bot committed
154
  writer->SetInput(filter->GetOutput());
155

156 157 158
  //  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.
159

160
  try
161
  {
162
    writer->Update();
163
  }
OTB Bot's avatar
STYLE  
OTB Bot committed
164
  catch (itk::ExceptionObject& excep)
165
  {
166 167
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
168
  }
OTB Bot's avatar
STYLE  
OTB Bot committed
169
  catch (...)
170
  {
171 172
    std::cout << "Unknown exception !" << std::endl;
    return EXIT_FAILURE;
173
  }
174 175

  // Pretty image creation for the printing
176 177 178 179 180 181 182 183 184 185 186
  typedef otb::Image<unsigned char, Dimension>                                                OutputPrettyImageType;
  typedef otb::VectorImage<unsigned char, Dimension>                                          OutputVectorPrettyImageType;
  typedef otb::ImageFileWriter<OutputVectorPrettyImageType>                                   WriterVectorPrettyType;
  typedef otb::ImageFileWriter<OutputPrettyImageType>                                         WriterPrettyType;
  typedef itk::RescaleIntensityImageFilter<OutputImageType, OutputPrettyImageType>            RescalerType;
  typedef otb::VectorRescaleIntensityImageFilter<InputImageType, OutputVectorPrettyImageType> VectorRescalerType;
  typedef otb::MultiChannelExtractROI<unsigned char, unsigned char>                           ChannelExtractorType;

  VectorRescalerType::Pointer     vectRescaler     = VectorRescalerType::New();
  ChannelExtractorType::Pointer   selecter         = ChannelExtractorType::New();
  WriterVectorPrettyType::Pointer vectPrettyWriter = WriterVectorPrettyType::New();
OTB Bot's avatar
STYLE  
OTB Bot committed
187 188

  OutputVectorPrettyImageType::PixelType minimum, maximum;
189 190 191 192 193 194
  minimum.SetSize(reader->GetOutput()->GetNumberOfComponentsPerPixel());
  maximum.SetSize(reader->GetOutput()->GetNumberOfComponentsPerPixel());
  minimum.Fill(0);
  maximum.Fill(255);
  vectRescaler->SetOutputMinimum(minimum);
  vectRescaler->SetOutputMaximum(maximum);
195
  vectRescaler->SetClampThreshold(0.01);
OTB Bot's avatar
STYLE  
OTB Bot committed
196
  vectRescaler->SetInput(reader->GetOutput());
197 198 199 200 201

  selecter->SetInput(vectRescaler->GetOutput());
  selecter->SetChannel(1);
  selecter->SetChannel(2);
  selecter->SetChannel(3);
202

OTB Bot's avatar
STYLE  
OTB Bot committed
203 204
  vectPrettyWriter->SetFileName(argv[3]);
  vectPrettyWriter->SetInput(selecter->GetOutput());
205

OTB Bot's avatar
STYLE  
OTB Bot committed
206
  typedef itk::ThresholdImageFilter<OutputImageType> ThresholderType;
207 208

  ThresholderType::Pointer thresholder = ThresholderType::New();
OTB Bot's avatar
STYLE  
OTB Bot committed
209 210 211
  thresholder->SetInput(filter->GetOutput());
  thresholder->SetOutsideValue(1.0);
  thresholder->ThresholdOutside(0.0, 1.0);
212 213
  thresholder->Update();

214
  RescalerType::Pointer     rescaler     = RescalerType::New();
215
  WriterPrettyType::Pointer prettyWriter = WriterPrettyType::New();
OTB Bot's avatar
STYLE  
OTB Bot committed
216
  rescaler->SetInput(thresholder->GetOutput());
217 218
  rescaler->SetOutputMinimum(0);
  rescaler->SetOutputMaximum(255);
OTB Bot's avatar
STYLE  
OTB Bot committed
219 220
  prettyWriter->SetFileName(argv[4]);
  prettyWriter->SetInput(rescaler->GetOutput());
221

222
  try
223
  {
224 225
    prettyWriter->Update();
    vectPrettyWriter->Update();
226
  }
OTB Bot's avatar
STYLE  
OTB Bot committed
227
  catch (itk::ExceptionObject& excep)
228
  {
229 230
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
231
  }
OTB Bot's avatar
STYLE  
OTB Bot committed
232
  catch (...)
233
  {
234 235
    std::cout << "Unknown exception !" << std::endl;
    return EXIT_FAILURE;
236
  }
237

238
  return EXIT_SUCCESS;
239

240
  // Let's now run this example using as input the image
241 242 243
  // \code{IndexVegetation.hd} (image kindly and free of charge given
  // by SISA and CNES) and $\gamma$=0.6 provided in the
  // directory \code{Examples/Data}.
244 245 246 247 248 249
  //
  //
  // \begin{figure} \center
  // \includegraphics[width=0.24\textwidth]{pretty_VegetationIndex.eps}
  // \includegraphics[width=0.24\textwidth]{pretty_ARVIMultiChannelRAndBAndNIRVegetationIndex.eps}
  // \itkcaption[ARVI Example]{ARVI result on the right with the left image in input.}
250
  // \label{fig:ARVIMultiChannelRAndBAndNIRIndexImageFilter}
251 252
  // \end{figure}
}