BayesianFusionImageFilter.cxx 12.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
/*
 * Copyright (C) 1999-2011 Insight Software Consortium
 * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
 *
 * 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.
 */
21

22

23 24 25 26 27

//  Software Guide : BeginCommandLineArgs
//  INPUTS: {multiSpect.tif} , {multiSpectInterp.tif}, {panchro.tif}
//  OUTPUTS: {BayesianFusion_0.9999.tif} , {pretty_BayesianFusion_0.9999.png} , {pretty_multiSpect_0.9999.png} , {pretty_multiSpectInterp_0.9999.png} , {pretty_panchro_0.9999.png}
//  0.9999
28
//  Software Guide : EndCommandLineArgs
29 30 31 32 33 34 35 36 37 38 39 40 41

//  Software Guide : BeginCommandLineArgs
//  INPUTS: {multiSpect.tif} , {multiSpectInterp.tif}, {panchro.tif}
//  OUTPUTS: {BayesianFusion_0.5.tif} , {pretty_BayesianFusion_0.5.png} , {pretty_multiSpect_0.5.png} , {pretty_multiSpectInterp_0.5.png} , {pretty_panchro_0.5.png}
//  0.5
//  Software Guide : EndCommandLineArgs

// Software Guide : BeginLatex
//
// \index{otb::BayesianFusionFilter}
// \index{otb::BayesianFusionFilter!header}
//
// The following example illustrates the use of the
42 43 44
// \doxygen{otb}{BayesianFusionFilter}. The Bayesian data fusion
// relies on the idea that variables of interest, denoted as vector $\mathbf{Z}$,
// cannot be directly observed. They are linked to the observable variables
45 46 47
// $\mathbf{Y}$ through the following error-like model.
//
// \begin{equation}
Jordi Inglada's avatar
Jordi Inglada committed
48
// \mathbf{Y} = g(\mathbf{Z}) + \mathbf{E}
49 50
// \end{equation}
//
51
// where g($\mathbf{Z}$) is a set of functionals and $\mathbf{E}$ is a
52 53
// vector of random errors that are stochastically independent from $\mathbf{Z}$.
// This algorithm uses elementary probability calculus, and several assumptions to compute
54
// the data fusion. For more explication see Fasbender, Radoux and Bogaert's
55 56 57 58 59 60
// publication \cite{JRadoux}.
// Three images are used :
// \begin{itemize}
// \item a panchromatic image,
// \item a multispectral image resampled at the panchromatic image spatial resolution,
// \item a multispectral image resampled at the panchromatic image spatial resolution,
61
// using, e.g. a cubic interpolator.
62 63 64 65
// \item a float : $\lambda$, the meaning of the weight to be given to the panchromatic
// image compared to the multispectral one.
// \end{itemize}
//
66
// Let's look at the minimal code required to use this algorithm. First, the following header
67
// defining the otb::BayesianFusionFilter class must be included.
68
// Software Guide : EndLatex
69 70 71 72 73 74 75 76 77 78 79 80 81

// Software Guide : BeginCodeSnippet
#include "otbBayesianFusionFilter.h"
// Software Guide : EndCodeSnippet

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

OTB Bot's avatar
OTB Bot committed
82
int main(int argc, char *argv[])
83
{
OTB Bot's avatar
OTB Bot committed
84 85
  if (argc < 10)
    {
86 87
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
OTB Bot's avatar
OTB Bot committed
88 89 90 91
    std::cerr << " inputMultiSpectralImage inputMultiSpectralInterpolatedImage "
              << "inputPanchromatiqueImage outputImage outputImagePrinted "
              << "msPrinted msiPrinted panchroPrinted lambda"
              << std::endl;
92
    return 1;
OTB Bot's avatar
OTB Bot committed
93
    }
94 95

  //  Software Guide : BeginLatex
96
  //
97 98 99 100
  //  The image types are now defined using pixel types and particular
  //  dimension. The panchromatic image is defined as an \doxygen{otb}{Image}
  //  and the multispectral one as \doxygen{otb}{VectorImage}.
  //
101
  //  Software Guide : EndLatex
102 103

  // Software Guide : BeginCodeSnippet
OTB Bot's avatar
OTB Bot committed
104 105 106 107
  typedef double InternalPixelType;
  const unsigned int Dimension = 2;
  typedef otb::Image<InternalPixelType, Dimension>       PanchroImageType;
  typedef otb::VectorImage<InternalPixelType, Dimension> MultiSpecImageType;
108 109
  // Software Guide : EndCodeSnippet

OTB Bot's avatar
OTB Bot committed
110 111
  typedef double                                       OutputPixelType;
  typedef otb::VectorImage<OutputPixelType, Dimension> OutputImageType;
112 113 114

  // We instantiate reader and writer types
  //
OTB Bot's avatar
OTB Bot committed
115 116 117
  typedef  otb::ImageFileReader<MultiSpecImageType> ReaderVectorType;
  typedef  otb::ImageFileReader<PanchroImageType>   ReaderType;
  typedef  otb::ImageFileWriter<OutputImageType>    WriterType;
118 119 120

  ReaderVectorType::Pointer multiSpectReader       = ReaderVectorType::New();
  ReaderVectorType::Pointer multiSpectInterpReader = ReaderVectorType::New();
OTB Bot's avatar
OTB Bot committed
121 122
  ReaderType::Pointer       panchroReader                = ReaderType::New();
  WriterType::Pointer       writer                       = WriterType::New();
123

OTB Bot's avatar
OTB Bot committed
124 125 126 127
  multiSpectReader->SetFileName(argv[1]);
  multiSpectInterpReader->SetFileName(argv[2]);
  panchroReader->SetFileName(argv[3]);
  writer->SetFileName(argv[4]);
128 129

  //  Software Guide : BeginLatex
130
  //
131 132 133
  //  The Bayesian data fusion filter type is instantiated using the images types as
  //  a template parameters.
  //
134
  //  Software Guide : EndLatex
135 136

  // Software Guide : BeginCodeSnippet
OTB Bot's avatar
OTB Bot committed
137
  typedef otb::BayesianFusionFilter<MultiSpecImageType,
OTB Bot's avatar
OTB Bot committed
138 139 140
      MultiSpecImageType,
      PanchroImageType,
      OutputImageType>
OTB Bot's avatar
OTB Bot committed
141
  BayesianFusionFilterType;
142 143 144
  // Software Guide : EndCodeSnippet

  //  Software Guide : BeginLatex
145
  //
146 147 148
  //  Next the filter is created by invoking the \code{New()} method and
  //  assigning the result to a \doxygen{itk}{SmartPointer}.
  //
149
  //  Software Guide : EndLatex
150 151

  // Software Guide : BeginCodeSnippet
OTB Bot's avatar
OTB Bot committed
152 153
  BayesianFusionFilterType::Pointer bayesianFilter =
    BayesianFusionFilterType::New();
154 155 156
  // Software Guide : EndCodeSnippet

  //  Software Guide : BeginLatex
157
  //
158
  //  Now the multi spectral image, the interpolated multi spectral image and
159
  //  the panchromatic image are given as inputs to the filter.
160
  //
161
  //  Software Guide : EndLatex
162 163

  // Software Guide : BeginCodeSnippet
OTB Bot's avatar
OTB Bot committed
164 165 166
  bayesianFilter->SetMultiSpect(multiSpectReader->GetOutput());
  bayesianFilter->SetMultiSpectInterp(multiSpectInterpReader->GetOutput());
  bayesianFilter->SetPanchro(panchroReader->GetOutput());
167

OTB Bot's avatar
OTB Bot committed
168
  writer->SetInput(bayesianFilter->GetOutput());
169 170 171 172
  // Software Guide : EndCodeSnippet

  //  Software Guide : BeginLatex
  //  The BayesianFusionFilter requires defining one parameter : $\lambda$.
173 174
  //  The $\lambda$ parameter can be used to tune the fusion toward either a high color
  //  consistency or sharp details. Typical $\lambda$ value range in  $[0.5, 1[$,  where higher
175 176
  //  values yield sharper details. by default $\lambda$ is set at 0.9999.
  //
177
  //  Software Guide : EndLatex
178 179

  // Software Guide : BeginCodeSnippet
OTB Bot's avatar
OTB Bot committed
180
  bayesianFilter->SetLambda(atof(argv[9]));
181 182 183
  // Software Guide : EndCodeSnippet

  //  Software Guide : BeginLatex
184
  //
185 186 187 188
  //  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.
  //
189
  //  Software Guide : EndLatex
190 191 192

  // Software Guide : BeginCodeSnippet
  try
OTB Bot's avatar
OTB Bot committed
193
    {
194
    writer->Update();
OTB Bot's avatar
OTB Bot committed
195 196 197
    }
  catch (itk::ExceptionObject& excep)
    {
198 199
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
OTB Bot's avatar
OTB Bot committed
200
    }
201
  // Software Guide : EndCodeSnippet
202

203
  // Create an 3 band images for the software guide
OTB Bot's avatar
OTB Bot committed
204 205 206
  typedef unsigned char                                 OutputPixelType2;
  typedef otb::VectorImage<OutputPixelType2, Dimension> OutputVectorImageType;
  typedef otb::ImageFileWriter<OutputVectorImageType>   VectorWriterType;
207
  typedef otb::VectorRescaleIntensityImageFilter<MultiSpecImageType,
OTB Bot's avatar
OTB Bot committed
208
      OutputVectorImageType>
OTB Bot's avatar
OTB Bot committed
209
  VectorRescalerType;
210
  typedef otb::VectorRescaleIntensityImageFilter<OutputImageType,
OTB Bot's avatar
OTB Bot committed
211
      OutputVectorImageType>
OTB Bot's avatar
OTB Bot committed
212 213
  VectorRescalerBayesianType;
  typedef otb::ImageToVectorImageCastFilter<PanchroImageType,
OTB Bot's avatar
OTB Bot committed
214
      MultiSpecImageType> CasterType;
OTB Bot's avatar
OTB Bot committed
215
  typedef otb::MultiChannelExtractROI<OutputPixelType2,
OTB Bot's avatar
OTB Bot committed
216 217
      OutputPixelType2>
  ChannelExtractorType;
218

219 220
  multiSpectReader->GenerateOutputInformation();
  multiSpectInterpReader->GenerateOutputInformation();
221

222
  CasterType::Pointer cast = CasterType::New();
223
  cast->SetInput(panchroReader->GetOutput());
224

OTB Bot's avatar
OTB Bot committed
225
  OutputVectorImageType::PixelType minimum, maximum;
226 227 228 229
  minimum.SetSize(multiSpectReader->GetOutput()->GetNumberOfComponentsPerPixel());
  maximum.SetSize(multiSpectReader->GetOutput()->GetNumberOfComponentsPerPixel());
  minimum.Fill(0);
  maximum.Fill(255);
230

231 232 233
  VectorRescalerType::Pointer         vrms  = VectorRescalerType::New();
  VectorRescalerType::Pointer         vrmsi = VectorRescalerType::New();
  VectorRescalerBayesianType::Pointer vrb   = VectorRescalerBayesianType::New();
234

235 236 237 238
  vrms->SetInput(multiSpectReader->GetOutput());
  vrms->SetOutputMinimum(minimum);
  vrms->SetOutputMaximum(maximum);
  vrms->SetClampThreshold(0.01);
239

240 241 242 243
  vrmsi->SetInput(multiSpectInterpReader->GetOutput());
  vrmsi->SetOutputMinimum(minimum);
  vrmsi->SetOutputMaximum(maximum);
  vrmsi->SetClampThreshold(0.01);
244

245 246 247 248
  vrb->SetInput(bayesianFilter->GetOutput());
  vrb->SetOutputMinimum(minimum);
  vrb->SetOutputMaximum(maximum);
  vrb->SetClampThreshold(0.01);
249

250 251 252 253 254 255 256 257 258
  VectorRescalerType::Pointer rp = VectorRescalerType::New();
  rp->SetInput(cast->GetOutput());
  minimum.SetSize(1);
  maximum.SetSize(1);
  minimum.Fill(0);
  maximum.Fill(255);
  rp->SetOutputMinimum(minimum);
  rp->SetOutputMaximum(maximum);
  rp->SetClampThreshold(0.01);
259

260 261 262 263 264
  ChannelExtractorType::Pointer selecterms = ChannelExtractorType::New();
  ChannelExtractorType::Pointer selectermsi = ChannelExtractorType::New();
  ChannelExtractorType::Pointer selecterf = ChannelExtractorType::New();

  selecterms->SetInput(vrms->GetOutput());
265
// selecterms->SetExtractionRegion(multiSpectReader->GetOutput()->GetLargestPossibleRegion());
266 267 268
  selecterms->SetChannel(2);
  selecterms->SetChannel(3);
  selecterms->SetChannel(4);
269

270
  selectermsi->SetInput(vrmsi->GetOutput());
271
// selectermsi->SetExtractionRegion(multiSpectInterpReader->GetOutput()->GetLargestPossibleRegion());
272 273 274
  selectermsi->SetChannel(2);
  selectermsi->SetChannel(3);
  selectermsi->SetChannel(4);
275

276 277 278 279 280
  selecterf->SetInput(vrb->GetOutput());
  //selecterf->SetExtractionRegion(bayesianFilter->GetOutput()->GetLargestPossibleRegion());
  selecterf->SetChannel(2);
  selecterf->SetChannel(3);
  selecterf->SetChannel(4);
281

282 283 284 285 286 287 288 289 290 291 292 293 294
  VectorWriterType::Pointer vectWriterms  = VectorWriterType::New();
  VectorWriterType::Pointer vectWritermsi = VectorWriterType::New();
  VectorWriterType::Pointer vectWriterf   = VectorWriterType::New();
  VectorWriterType::Pointer vectWriterp   = VectorWriterType::New();

  vectWriterf->SetFileName(argv[5]);
  vectWriterf->SetInput(selecterf->GetOutput());
  vectWriterms->SetFileName(argv[6]);
  vectWriterms->SetInput(selecterms->GetOutput());
  vectWritermsi->SetFileName(argv[7]);
  vectWritermsi->SetInput(selectermsi->GetOutput());
  vectWriterp->SetFileName(argv[8]);
  vectWriterp->SetInput(rp->GetOutput());
295

296
  try
OTB Bot's avatar
OTB Bot committed
297
    {
298 299 300 301
    vectWriterms->Update();
    vectWritermsi->Update();
    vectWriterf->Update();
    vectWriterp->Update();
OTB Bot's avatar
OTB Bot committed
302 303 304
    }
  catch (itk::ExceptionObject& excep)
    {
305 306
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
OTB Bot's avatar
OTB Bot committed
307 308 309
    }
  catch (...)
    {
310 311
    std::cout << "Unknown exception !" << std::endl;
    return EXIT_FAILURE;
OTB Bot's avatar
OTB Bot committed
312
    }
313 314 315 316 317

  //  Software Guide : BeginLatex
  //
  //  Let's now run this example using as input the images
  //  \code{multiSpect.tif} , \code{multiSpectInterp.tif} and \code{panchro.tif}
318 319 320
  //  provided in the directory \code{Examples/Data}. The results
  //  obtained for 2 different values for $\lambda$ are shown in figure
  //  \ref{fig:BayesianImageFusionFilterInput}.
321 322 323 324 325 326 327
  //
  //
  // \begin{figure} \center
  // \includegraphics[width=0.24\textwidth]{pretty_multiSpect_0.5.eps}
  // \includegraphics[width=0.24\textwidth]{pretty_multiSpectInterp_0.5.eps}
  // \includegraphics[width=0.24\textwidth]{pretty_panchro_0.5.eps}
  // \itkcaption[Bayesian Data Fusion Example inputs]{Input
328
  // images used for this example (\copyright European Space Imaging).}
329 330
  // \label{fig:BayesianImageFusionFilterInput}
  // \end{figure}
331

332 333 334 335 336 337 338 339 340
  // \begin{figure} \center
  // \includegraphics[width=0.24\textwidth]{pretty_BayesianFusion_0.5.eps}
  // \includegraphics[width=0.24\textwidth]{pretty_BayesianFusion_0.9999.eps}
  // \itkcaption[Bayesian Data Fusion results]{Fusion results
  // for the Bayesian Data Fusion filter for $\lambda = 0.5$ on the left and $\lambda = 0.9999$ on the right.}
  // \label{fig:BayesianImageFusionFilterOutput}
  // \end{figure}
  //
  //
341
  //  Software Guide : EndLatex
342

343
  return EXIT_SUCCESS;
344
}