NCCRegistrationFilterExample.cxx 7.75 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:
./NCCRegistrationFilterExample Input/StereoFixed.png \
                               Input/StereoMoving.png \
                               Output/deformationFieldOutput-horizontal.png \
                               Output/deformationFieldOutput-vertical.png \
                               Output/resampledOutput2.png \
                               5 \
                               1.0 \
                               2
*/
32 33


34
// This example demonstrates the use of the \doxygen{otb}{NCCRegistrationFilter}. This filter performs deformation estimation
35
// by optimising a PDE based on the normalized correlation coefficient. It uses the finite difference solver hierarchy.
36 37
//
// The first step toward the use of these filters is to include the proper header files.
38

39
#include "otbImageFileWriter.h"
40
#include "otbImageFileReader.h"
41 42

#include "otbNCCRegistrationFilter.h"
43
#include "itkRecursiveGaussianImageFilter.h"
44
#include "itkWarpImageFilter.h"
45

46
#include "otbImageOfVectorsToMonoChannelExtractROI.h"
47
#include "itkUnaryFunctorImageFilter.h"
48 49 50
#include "itkRescaleIntensityImageFilter.h"
#include "itkCastImageFilter.h"

51 52
#include <iostream>

OTB Bot's avatar
STYLE  
OTB Bot committed
53
int main(int argc, char** argv)
54
{
55

OTB Bot's avatar
STYLE  
OTB Bot committed
56
  if (argc != 9)
57
  {
OTB Bot's avatar
STYLE  
OTB Bot committed
58
    std::cerr << "Usage: " << argv[0];
59
    std::cerr << " fixedFileName movingFileName fieldOutNameHorizontal fieldOutNameVertical imageOutName ";
OTB Bot's avatar
STYLE  
OTB Bot committed
60
    std::cerr << "explorationSize bluringSigma nbIterations ";
61

62
    return EXIT_FAILURE;
63
  }
64

65 66
  const unsigned int ImageDimension = 2;

OTB Bot's avatar
STYLE  
OTB Bot committed
67
  typedef double                              PixelType;
68
  typedef itk::Vector<double, ImageDimension> DisplacementPixelType;
69

OTB Bot's avatar
STYLE  
OTB Bot committed
70 71
  typedef unsigned char                               OutputPixelType;
  typedef otb::Image<OutputPixelType, ImageDimension> OutputImageType;
72

73
  // Several type of \doxygen{otb}{Image} are required to represent the reference image (fixed)
74
  // the image we want to register (moving) and the deformation field.
75

76 77 78 79
  // Allocate Images
  typedef otb::Image<PixelType, ImageDimension>             MovingImageType;
  typedef otb::Image<PixelType, ImageDimension>             FixedImageType;
  typedef otb::Image<DisplacementPixelType, ImageDimension> DisplacementFieldType;
80

OTB Bot's avatar
STYLE  
OTB Bot committed
81
  typedef otb::ImageFileReader<FixedImageType> FixedReaderType;
82
  FixedReaderType::Pointer                     fReader = FixedReaderType::New();
83
  fReader->SetFileName(argv[1]);
84

OTB Bot's avatar
STYLE  
OTB Bot committed
85
  typedef otb::ImageFileReader<MovingImageType> MovingReaderType;
86
  MovingReaderType::Pointer                     mReader = MovingReaderType::New();
87
  mReader->SetFileName(argv[2]);
88

89 90 91
  // To make the correlation estimation more robust, the first
  // required step is to blur the input images. This is done using the
  // \doxygen{itk}{RecursiveGaussianImageFilter}:
92

93 94
  // Blur input images
  typedef itk::RecursiveGaussianImageFilter<FixedImageType, FixedImageType> FixedBlurType;
95 96

  FixedBlurType::Pointer fBlur = FixedBlurType::New();
OTB Bot's avatar
STYLE  
OTB Bot committed
97
  fBlur->SetInput(fReader->GetOutput());
98
  fBlur->SetSigma(std::stof(argv[7]));
99

100
  typedef itk::RecursiveGaussianImageFilter<MovingImageType, MovingImageType> MovingBlurType;
101 102

  MovingBlurType::Pointer mBlur = MovingBlurType::New();
OTB Bot's avatar
STYLE  
OTB Bot committed
103
  mBlur->SetInput(mReader->GetOutput());
104
  mBlur->SetSigma(std::stof(argv[7]));
105

106
  // Now, we need to instantiate the NCCRegistrationFilter which is going to perform the registration:
107

108 109
  // Create the filter
  typedef otb::NCCRegistrationFilter<FixedImageType, MovingImageType, DisplacementFieldType> RegistrationFilterType;
110 111 112

  RegistrationFilterType::Pointer registrator = RegistrationFilterType::New();

OTB Bot's avatar
STYLE  
OTB Bot committed
113 114
  registrator->SetMovingImage(mBlur->GetOutput());
  registrator->SetFixedImage(fBlur->GetOutput());
115

116
  // Some parameters need to be specified to the NCCRegistrationFilter:
117 118
  // \begin{itemize}
  // \item The area where the search is performed. This area is defined by its radius:
119

120 121 122 123
  typedef RegistrationFilterType::RadiusType RadiusType;

  RadiusType radius;

124 125
  radius[0] = std::stoi(argv[6]);
  radius[1] = std::stoi(argv[6]);
126

OTB Bot's avatar
STYLE  
OTB Bot committed
127
  registrator->SetNCCRadius(radius);
128

129
  std::cout << "NCC radius " << registrator->GetNCCRadius() << std::endl;
130

131
  // \item The number of iterations for the PDE resolution:
132

133
  registrator->SetNumberOfIterations(std::stoi(argv[8]));
134
  // registrator->GetDisplacementField();
135

136
  // \end{itemize}
137 138 139
  // The execution of the NCCRegistrationFilter will be triggered by
  // the \code{Update()} call on the writer at the end of the
  // pipeline. Make sure to use a
140
  // \doxygen{otb}{ImageFileWriter} if you want to benefit
141
  // from the streaming features.
142

143 144
  typedef otb::ImageOfVectorsToMonoChannelExtractROI<DisplacementFieldType, MovingImageType> ChannelExtractionFilterType;
  ChannelExtractionFilterType::Pointer                                                       channelExtractor = ChannelExtractionFilterType::New();
145 146 147

  channelExtractor->SetInput(registrator->GetOutput());
  channelExtractor->SetChannel(1);
148

149 150
  typedef itk::RescaleIntensityImageFilter<MovingImageType, OutputImageType> RescalerType;
  RescalerType::Pointer                                                      fieldRescaler = RescalerType::New();
151

152 153 154
  fieldRescaler->SetInput(channelExtractor->GetOutput());
  fieldRescaler->SetOutputMaximum(255);
  fieldRescaler->SetOutputMinimum(0);
155

156
  typedef otb::ImageFileWriter<OutputImageType> DFWriterType;
157 158 159

  DFWriterType::Pointer dfWriter = DFWriterType::New();
  dfWriter->SetFileName(argv[3]);
OTB Bot's avatar
STYLE  
OTB Bot committed
160
  dfWriter->SetInput(fieldRescaler->GetOutput());
161 162
  dfWriter->Update();

163 164 165
  channelExtractor->SetChannel(2);
  dfWriter->SetFileName(argv[4]);
  dfWriter->Update();
166

167 168
  typedef itk::WarpImageFilter<MovingImageType, MovingImageType, DisplacementFieldType> WarperType;
  WarperType::Pointer                                                                   warper = WarperType::New();
169

170 171
  MovingImageType::PixelType padValue = 4.0;

OTB Bot's avatar
STYLE  
OTB Bot committed
172
  warper->SetInput(mReader->GetOutput());
173
  warper->SetDisplacementField(registrator->GetOutput());
OTB Bot's avatar
STYLE  
OTB Bot committed
174
  warper->SetEdgePaddingValue(padValue);
175

OTB Bot's avatar
STYLE  
OTB Bot committed
176
  typedef itk::CastImageFilter<MovingImageType, OutputImageType> CastFilterType;
177
  CastFilterType::Pointer                                        caster = CastFilterType::New();
OTB Bot's avatar
STYLE  
OTB Bot committed
178
  caster->SetInput(warper->GetOutput());
179

180
  typedef otb::ImageFileWriter<OutputImageType> WriterType;
181 182

  WriterType::Pointer writer = WriterType::New();
183
  writer->SetFileName(argv[5]);
OTB Bot's avatar
STYLE  
OTB Bot committed
184
  writer->SetInput(caster->GetOutput());
185
  writer->Update();
186

Jordi Inglada's avatar
Jordi Inglada committed
187 188
  // Figure~\ref{fig:NCCRegistrationFilterOUTPUT} shows the result of
  // applying the disparity map estimation.
189
  //
190
  // \begin{figure}
191
  // \center
192 193
  // \includegraphics[width=0.40\textwidth]{StereoFixed.eps}
  // \includegraphics[width=0.40\textwidth]{StereoMoving.eps}
194 195
  // \includegraphics[width=0.40\textwidth]{deformationFieldOutput-horizontal.eps}
  // \includegraphics[width=0.40\textwidth]{deformationFieldOutput-vertical.eps}
196
  // \itkcaption[Displacement field and resampling from NCC registration]{From left
197
  // to right and top to bottom: fixed input image, moving image with a low stereo angle,
198
  // estimated deformation field in the horizontal direction, estimated deformation field in the vertical direction.}
199 200
  // \label{fig:NCCRegistrationFilterOUTPUT}
  // \end{figure}
201

202 203
  return EXIT_SUCCESS;
}