StereoReconstructionExample.cxx 17.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/*
 * 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.
 */
20 21 22 23 24



//  Software Guide : BeginCommandLineArgs
//    INPUTS: {sensor_stereo_left.tif}, {sensor_stereo_right.tif}
25
//    OUTPUTS: {elevationOutput.tif},{elevationOutputPrintable.png}
OTB Bot's avatar
STYLE  
OTB Bot committed
26
//    140
27 28 29 30
//  Software Guide : EndCommandLineArgs

// Software Guide : BeginLatex
//
OTB Bot's avatar
STYLE  
OTB Bot committed
31 32 33
// This example demonstrates the use of the stereo reconstruction chain
// from an image pair. The images are assumed to come from the same sensor
// but with different positions. The approach presented here has the
34
// following steps:
35
// \begin{itemize}
36 37
// \item Epipolar resampling of the image pair
// \item Dense disparity map estimation
38
// \item Projection of the disparities on an existing Digital Elevation Model (DEM)
39 40 41 42 43
// \end{itemize}
// It is important to note that this method requires the sensor models with
// a pose estimate for each image.


44 45 46 47
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
48
#include "otbStereorectificationDisplacementFieldSource.h"
49
#include "otbStreamingWarpImageFilter.h"
50
#include "otbBandMathImageFilter.h"
51 52 53 54 55
#include "otbSubPixelDisparityImageFilter.h"
#include "otbDisparityMapMedianFilter.h"
#include "otbDisparityMapToDEMFilter.h"

#include "otbImageFileReader.h"
56
#include "otbImageFileWriter.h"
57
#include "otbBCOInterpolateImageFunction.h"
58
#include "itkUnaryFunctorImageFilter.h"
59 60 61
#include "itkVectorCastImageFilter.h"
#include "otbImageList.h"
#include "otbImageListToVectorImageFilter.h"
62
#include "itkRescaleIntensityImageFilter.h"
63
#include "otbDEMHandler.h"
64 65 66 67 68
// Software Guide : EndCodeSnippet


int main(int argc, char* argv[])
{
69
  if (argc != 6)
70 71
    {
    std::cerr << "Usage: " << argv[0];
72
    std::cerr << " sensorImage1 sensorImage2 outputDEM outputDEMPNG";
73 74 75
    std::cerr << "averageElevation  " << std::endl;
    return EXIT_FAILURE;
    }
76

77 78
  typedef otb::Image<float,2>                 FloatImageType;
  typedef otb::VectorImage<float,2>           FloatVectorImageType;
79

80 81
  typedef otb::ImageFileReader
    <FloatImageType>                          ImageReaderType;
82

83
  typedef otb::ImageFileWriter
84
    <FloatImageType>                          WriterType;
85

86 87 88 89 90 91
  typedef unsigned char                       OutputPixelType;
  typedef otb::Image<OutputPixelType, 2>      OutputImageType;

  typedef itk::RescaleIntensityImageFilter<FloatImageType,
      OutputImageType> RescalerType;

92
  typedef otb::ImageFileWriter
93
    <OutputImageType>                         OutputWriterType;
94 95 96
// Software Guide : BeginLatex
// This example demonstrates the use of the following filters :
// \begin{itemize}
97
// \item \doxygen{otb}{StereorectificationDisplacementFieldSource}
98 99 100 101 102
// \item \doxygen{otb}{StreamingWarpImageFilter}
// \item \doxygen{otb}{PixelWiseBlockMatchingImageFilter}
// \item \doxygen{otb}{otbSubPixelDisparityImageFilter}
// \item \doxygen{otb}{otbDisparityMapMedianFilter}
// \item \doxygen{otb}{DisparityMapToDEMFilter}
103
// \end{itemize}
104 105 106
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
107 108
  typedef otb::StereorectificationDisplacementFieldSource
    <FloatImageType,FloatVectorImageType>     DisplacementFieldSourceType;
109

110 111
  typedef itk::Vector<double,2>               DisplacementType;
  typedef otb::Image<DisplacementType>         DisplacementFieldType;
112

113 114
  typedef itk::VectorCastImageFilter
    <FloatVectorImageType,
115
     DisplacementFieldType>                    DisplacementFieldCastFilterType;
116

117 118 119
  typedef otb::StreamingWarpImageFilter
    <FloatImageType,
     FloatImageType,
120
     DisplacementFieldType>                    WarpFilterType;
121

122 123
  typedef otb::BCOInterpolateImageFunction
    <FloatImageType>                          BCOInterpolationType;
124

125 126
  typedef otb::Functor::NCCBlockMatching
    <FloatImageType,FloatImageType>           NCCBlockMatchingFunctorType;
127

128 129 130 131 132
  typedef otb::PixelWiseBlockMatchingImageFilter
    <FloatImageType,
     FloatImageType,
     FloatImageType,
     FloatImageType,
133
     NCCBlockMatchingFunctorType>             NCCBlockMatchingFilterType;
134

135 136
  typedef otb::BandMathImageFilter
    <FloatImageType>                          BandMathFilterType;
137

138 139 140 141 142
  typedef otb::SubPixelDisparityImageFilter
    <FloatImageType,
     FloatImageType,
     FloatImageType,
     FloatImageType,
143
     NCCBlockMatchingFunctorType>             NCCSubPixelDisparityFilterType;
144

145 146 147
  typedef otb::DisparityMapMedianFilter
    <FloatImageType,
     FloatImageType,
148
     FloatImageType>                          MedianFilterType;
149

150 151 152 153 154 155 156
  typedef otb::DisparityMapToDEMFilter
    <FloatImageType,
     FloatImageType,
     FloatImageType,
     FloatVectorImageType,
     FloatImageType>                          DisparityToElevationFilterType;
// Software Guide : EndCodeSnippet
157

158
  double avgElevation = atof(argv[5]);
159
  otb::DEMHandler::Instance()->SetDefaultHeightAboveEllipsoid(avgElevation);
160

161 162
  ImageReaderType::Pointer leftReader =  ImageReaderType::New();
  ImageReaderType::Pointer rightReader =  ImageReaderType::New();
163

164 165
  leftReader->SetFileName(argv[1]);
  rightReader->SetFileName(argv[2]);
166

167 168
// Software Guide : BeginLatex
// The image pair is supposed to be in sensor geometry. From two images covering
OTB Bot's avatar
STYLE  
OTB Bot committed
169
// nearly the same area, one can estimate a common epipolar geometry. In this geometry,
170
// an altitude variation corresponds to an horizontal shift between the two images.
171
// The filter \doxygen{otb}{StereorectificationDisplacementFieldSource} computes the
172
// deformation grids for each image.
OTB Bot's avatar
STYLE  
OTB Bot committed
173
//
174 175 176 177 178 179
// These grids are sampled in epipolar geometry. They have two bands, containing
// the position offset (in physical space units) between the current epipolar
// point and the corresponding sensor point in horizontal and vertical
// direction. They can be computed at a lower resolution than sensor
// resolution. The application \code{StereoRectificationGridGenerator} also
// provides a simple tool to generate the epipolar grids for your image pair.
180
// Software Guide : EndLatex
181

182
// Software Guide : BeginCodeSnippet
183 184 185 186 187 188 189 190
  DisplacementFieldSourceType::Pointer m_DisplacementFieldSource = DisplacementFieldSourceType::New();
  m_DisplacementFieldSource->SetLeftImage(leftReader->GetOutput());
  m_DisplacementFieldSource->SetRightImage(rightReader->GetOutput());
  m_DisplacementFieldSource->SetGridStep(4);
  m_DisplacementFieldSource->SetScale(1.0);
  //m_DisplacementFieldSource->SetAverageElevation(avgElevation);

  m_DisplacementFieldSource->Update();
191
// Software Guide : EndCodeSnippet
192

193
// Software Guide : BeginLatex
OTB Bot's avatar
STYLE  
OTB Bot committed
194 195
// Then, the sensor images can be resampled in epipolar geometry, using the
// \doxygen{otb}{StreamingWarpImageFilter}. The application
196
// \code{GridBasedImageResampling} also gives an easy access to this filter. The user
197 198
// can choose the epipolar region to resample, as well as the resampling step
// and the interpolator.
199 200 201 202 203 204
//
// Note that the epipolar image size can be retrieved from the stereo rectification grid
// filter.
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
205 206 207
  FloatImageType::SpacingType epipolarSpacing;
  epipolarSpacing[0] = 1.0;
  epipolarSpacing[1] = 1.0;
208

209
  FloatImageType::SizeType epipolarSize;
210
  epipolarSize = m_DisplacementFieldSource->GetRectifiedImageSize();
211

212 213 214
  FloatImageType::PointType epipolarOrigin;
  epipolarOrigin[0] = 0.0;
  epipolarOrigin[1] = 0.0;
215

216
  FloatImageType::PixelType defaultValue = 0;
217
// Software Guide : EndCodeSnippet
218

219
// Software Guide : BeginLatex
220
// The deformation grids are casted into deformation fields, then the left
221
// and right sensor images are resampled.
222
// Software Guide : EndLatex
223

224
// Software Guide : BeginCodeSnippet
225 226 227
  DisplacementFieldCastFilterType::Pointer m_LeftDisplacementFieldCaster = DisplacementFieldCastFilterType::New();
  m_LeftDisplacementFieldCaster->SetInput(m_DisplacementFieldSource->GetLeftDisplacementFieldOutput());
  m_LeftDisplacementFieldCaster->GetOutput()->UpdateOutputInformation();
228

229 230
  BCOInterpolationType::Pointer leftInterpolator = BCOInterpolationType::New();
  leftInterpolator->SetRadius(2);
231

232 233
  WarpFilterType::Pointer m_LeftWarpImageFilter = WarpFilterType::New();
  m_LeftWarpImageFilter->SetInput(leftReader->GetOutput());
234
  m_LeftWarpImageFilter->SetDisplacementField(m_LeftDisplacementFieldCaster->GetOutput());
235 236 237 238 239
  m_LeftWarpImageFilter->SetInterpolator(leftInterpolator);
  m_LeftWarpImageFilter->SetOutputSize(epipolarSize);
  m_LeftWarpImageFilter->SetOutputSpacing(epipolarSpacing);
  m_LeftWarpImageFilter->SetOutputOrigin(epipolarOrigin);
  m_LeftWarpImageFilter->SetEdgePaddingValue(defaultValue);
240

241 242 243
  DisplacementFieldCastFilterType::Pointer m_RightDisplacementFieldCaster = DisplacementFieldCastFilterType::New();
  m_RightDisplacementFieldCaster->SetInput(m_DisplacementFieldSource->GetRightDisplacementFieldOutput());
  m_RightDisplacementFieldCaster->GetOutput()->UpdateOutputInformation();
244

245 246
  BCOInterpolationType::Pointer rightInterpolator = BCOInterpolationType::New();
  rightInterpolator->SetRadius(2);
247

248 249
  WarpFilterType::Pointer m_RightWarpImageFilter = WarpFilterType::New();
  m_RightWarpImageFilter->SetInput(rightReader->GetOutput());
250
  m_RightWarpImageFilter->SetDisplacementField(m_RightDisplacementFieldCaster->GetOutput());
251 252 253 254 255
  m_RightWarpImageFilter->SetInterpolator(rightInterpolator);
  m_RightWarpImageFilter->SetOutputSize(epipolarSize);
  m_RightWarpImageFilter->SetOutputSpacing(epipolarSpacing);
  m_RightWarpImageFilter->SetOutputOrigin(epipolarOrigin);
  m_RightWarpImageFilter->SetEdgePaddingValue(defaultValue);
256
// Software Guide : EndCodeSnippet
257

258
// Software Guide : BeginLatex
OTB Bot's avatar
STYLE  
OTB Bot committed
259
// Since the resampling produces black regions around the image, it is useless
260
// to estimate disparities on these \textit{no-data} regions. We use a \doxygen{otb}{BandMathImageFilter}
261 262 263 264
// to produce a mask on left and right epipolar images.
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
265 266
  BandMathFilterType::Pointer m_LBandMathFilter = BandMathFilterType::New();
  m_LBandMathFilter->SetNthInput(0,m_LeftWarpImageFilter->GetOutput(),"inleft");
267 268 269
  #ifdef OTB_MUPARSER_HAS_CXX_LOGICAL_OPERATORS
  std::string leftExpr = "inleft != 0 ? 255 : 0";
  #else
270
  std::string leftExpr = "if(inleft != 0,255,0)";
271
  #endif
272

273
  m_LBandMathFilter->SetExpression(leftExpr);
274

275 276
  BandMathFilterType::Pointer m_RBandMathFilter = BandMathFilterType::New();
  m_RBandMathFilter->SetNthInput(0,m_RightWarpImageFilter->GetOutput(),"inright");
277 278 279 280

  #ifdef OTB_MUPARSER_HAS_CXX_LOGICAL_OPERATORS
  std::string rightExpr = "inright != 0 ? 255 : 0";
  #else
281
  std::string rightExpr = "if(inright != 0,255,0)";
282 283
  #endif

284
  m_RBandMathFilter->SetExpression(rightExpr);
285 286 287
// Software Guide : EndCodeSnippet

// Software Guide : BeginLatex
OTB Bot's avatar
STYLE  
OTB Bot committed
288 289 290 291 292 293
// Once the two sensor images have been resampled in epipolar geometry, the
// disparity map can be computed. The approach presented here is a 2D matching
// based on a pixel-wise metric optimization. This approach doesn't give the best
// results compared to global optimization methods, but it is suitable for
// streaming and threading on large images.
//
294
// The major filter used for this step is \doxygen{otb}{PixelWiseBlockMatchingImageFilter}.
OTB Bot's avatar
STYLE  
OTB Bot committed
295
// The metric is computed on a window centered around the tested epipolar position.
296
// It performs a pixel-to-pixel matching between the two epipolar images. The output disparities
OTB Bot's avatar
STYLE  
OTB Bot committed
297
// are given as index offset from left to right position. The following features are available
298 299 300 301 302 303
// in this filter:
// \begin{itemize}
// \item Available metrics : SSD, NCC and $L^{p}$ pseudo norm (computed on a square window)
// \item Rectangular disparity exploration area.
// \item Input masks for left and right images (optional).
// \item Output metric values (optional).
OTB Bot's avatar
STYLE  
OTB Bot committed
304
// \item Possibility to use input disparity estimate (as a uniform value or a full map) and an
305 306 307
// exploration radius around these values to reduce the size of the exploration area (optional).
// \end{itemize}
// Software Guide : EndLatex
308

309
// Software Guide : BeginCodeSnippet
310 311 312 313 314 315 316 317 318 319 320
  NCCBlockMatchingFilterType::Pointer m_NCCBlockMatcher = NCCBlockMatchingFilterType::New();
  m_NCCBlockMatcher->SetLeftInput(m_LeftWarpImageFilter->GetOutput());
  m_NCCBlockMatcher->SetRightInput(m_RightWarpImageFilter->GetOutput());
  m_NCCBlockMatcher->SetRadius(3);
  m_NCCBlockMatcher->SetMinimumHorizontalDisparity(-24);
  m_NCCBlockMatcher->SetMaximumHorizontalDisparity(0);
  m_NCCBlockMatcher->SetMinimumVerticalDisparity(0);
  m_NCCBlockMatcher->SetMaximumVerticalDisparity(0);
  m_NCCBlockMatcher->MinimizeOff();
  m_NCCBlockMatcher->SetLeftMaskInput(m_LBandMathFilter->GetOutput());
  m_NCCBlockMatcher->SetRightMaskInput(m_RBandMathFilter->GetOutput());
321
// Software Guide : EndCodeSnippet
322

323
// Software Guide : BeginLatex
324
// Some other filters have been added to enhance these \textit{pixel-to-pixel} disparities. The filter
OTB Bot's avatar
STYLE  
OTB Bot committed
325
// \doxygen{otb}{SubPixelDisparityImageFilter} can estimate the disparities with sub-pixel
326
// precision. Several interpolation methods can be used : parabolic fit, triangular fit, and
327 328 329 330
// dichotomy search.
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
331 332 333
  NCCSubPixelDisparityFilterType::Pointer m_NCCSubPixFilter = NCCSubPixelDisparityFilterType::New();
  m_NCCSubPixFilter->SetInputsFromBlockMatchingFilter(m_NCCBlockMatcher);
  m_NCCSubPixFilter->SetRefineMethod(NCCSubPixelDisparityFilterType::DICHOTOMY);
334 335 336 337 338 339 340
// Software Guide : EndCodeSnippet

// Software Guide : BeginLatex
// The filter \doxygen{otb}{DisparityMapMedianFilter} can be used to remove outliers. It has two
// parameters:
// \begin{itemize}
// \item The radius of the local neighborhood to compute the median
341
// \item An incoherence threshold to reject disparities whose distance from the local median
342 343 344
// is superior to the threshold.
// \end{itemize}
// Software Guide : EndLatex
345

346 347 348 349 350 351
// Software Guide : BeginCodeSnippet
  MedianFilterType::Pointer m_HMedianFilter = MedianFilterType::New();
  m_HMedianFilter->SetInput(m_NCCSubPixFilter->GetHorizontalDisparityOutput());
  m_HMedianFilter->SetRadius(2);
  m_HMedianFilter->SetIncoherenceThreshold(2.0);
  m_HMedianFilter->SetMaskInput(m_LBandMathFilter->GetOutput());
352

353 354 355 356 357 358 359 360
  MedianFilterType::Pointer m_VMedianFilter = MedianFilterType::New();
  m_VMedianFilter->SetInput(m_NCCSubPixFilter->GetVerticalDisparityOutput());
  m_VMedianFilter->SetRadius(2);
  m_VMedianFilter->SetIncoherenceThreshold(2.0);
  m_VMedianFilter->SetMaskInput(m_LBandMathFilter->GetOutput());
// Software Guide : EndCodeSnippet

// Software Guide : BeginLatex
361 362
// The application \code{PixelWiseBlockMatching} contains all these filters and
// provides a single interface to compute your disparity maps.
OTB Bot's avatar
STYLE  
OTB Bot committed
363
//
364 365
// The disparity map obtained with the previous step usually gives a good idea of
// the altitude profile. However, it is more useful to study altitude with a DEM (Digital
OTB Bot's avatar
STYLE  
OTB Bot committed
366 367
// Elevation Model) representation.
//
368
// The filter \doxygen{otb}{DisparityMapToDEMFilter} performs this last step. The behavior
369 370 371 372 373
// of this filter is to :
// \begin{itemize}
// \item Compute the DEM extent from the left sensor image envelope (spacing is set by the user)
// \item Compute the left and right rays corresponding to each valid disparity
// \item Compute the intersection with the \textit{mid-point} method
OTB Bot's avatar
STYLE  
OTB Bot committed
374
// \item If the 3D point falls inside a DEM cell and has a greater elevation than the
375 376
// current height, the cell height is updated
// \end{itemize}
OTB Bot's avatar
STYLE  
OTB Bot committed
377 378 379 380
// The rule of keeping the highest elevation makes sense for buildings seen from the side
// because the roof edges elevation has to be kept. However this rule is not suited for
// noisy disparities.
//
381 382 383 384 385 386 387 388 389
// The application \code{DisparityMapToElevationMap} also gives an example of use.
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
  DisparityToElevationFilterType::Pointer m_DispToElev = DisparityToElevationFilterType::New();
  m_DispToElev->SetHorizontalDisparityMapInput(m_HMedianFilter->GetOutput());
  m_DispToElev->SetVerticalDisparityMapInput(m_VMedianFilter->GetOutput());
  m_DispToElev->SetLeftInput(leftReader->GetOutput());
  m_DispToElev->SetRightInput(rightReader->GetOutput());
390 391
  m_DispToElev->SetLeftEpipolarGridInput(m_DisplacementFieldSource->GetLeftDisplacementFieldOutput());
  m_DispToElev->SetRightEpipolarGridInput(m_DisplacementFieldSource->GetRightDisplacementFieldOutput());
392 393
  m_DispToElev->SetElevationMin(avgElevation-10.0);
  m_DispToElev->SetElevationMax(avgElevation+80.0);
394 395
  m_DispToElev->SetDEMGridStep(2.5);
  m_DispToElev->SetDisparityMaskInput(m_LBandMathFilter->GetOutput());
396
  //m_DispToElev->SetAverageElevation(avgElevation);
397

398 399 400 401
  WriterType::Pointer m_DEMWriter = WriterType::New();
  m_DEMWriter->SetInput(m_DispToElev->GetOutput());
  m_DEMWriter->SetFileName(argv[3]);
  m_DEMWriter->Update();
402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427

  RescalerType::Pointer fieldRescaler = RescalerType::New();
  fieldRescaler->SetInput(m_DispToElev->GetOutput());
  fieldRescaler->SetOutputMaximum(255);
  fieldRescaler->SetOutputMinimum(0);

  OutputWriterType::Pointer fieldWriter = OutputWriterType::New();
  fieldWriter->SetInput(fieldRescaler->GetOutput());
  fieldWriter->SetFileName(argv[4]);
  fieldWriter->Update();
  // Software Guide : EndCodeSnippet

  // Software Guide : BeginLatex
  //
  // Figure~\ref{fig:STEREORECONSTRUCTIONOUTPUT} shows the result of applying
  // terrain reconstruction based using pixel-wise block matching, sub-pixel
  // interpolation and DEM estimation using a pair of Pleiades images over the
  // \textit{Stadium} in Toulouse, France.
  // \begin{figure}
  // \center
  // \includegraphics[width=0.4\textwidth]{elevationOutputPrintable.eps}
  // \itkcaption[From stereo pair to elevation]{DEM image estimated from the disparity.}
  // \label{fig:STEREORECONSTRUCTIONOUTPUT}
  // \end{figure}
  //
  // Software Guide : EndLatex
428

429 430
  return EXIT_SUCCESS;
}