otbImageToSIFTKeyPointSetFilterDistanceMap.cxx 12.6 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

#include <iostream>
#include <fstream>

24
#include "itkUnaryFunctorImageFilter.h"
25 26 27 28 29 30 31 32 33 34 35
#include "itkPointSet.h"
#include "itkVariableLengthVector.h"
#include "itkResampleImageFilter.h"
#include "itkDanielssonDistanceMapImageFilter.h"
#include "itkPointSetToImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"

#include "otbImageToSIFTKeyPointSetFilter.h"
#include "otbImage.h"
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
36
#include "itkAffineTransform.h"
37 38

typedef itk::VariableLengthVector<float> RealVectorType;
OTB Bot's avatar
STYLE  
OTB Bot committed
39 40 41
typedef itk::PointSet<RealVectorType, 2> PointSetType;
typedef otb::Image<float, 2>             ImageType;
typedef otb::Image<unsigned char, 2>     OutputImageType;
42

43
// PointSet iterator types
OTB Bot's avatar
STYLE  
OTB Bot committed
44 45
typedef PointSetType::PointsContainer    PointsContainerType;
typedef PointsContainerType::Iterator    PointsIteratorType;
46 47 48 49
typedef PointSetType::PointDataContainer PointDataContainerType;
typedef PointDataContainerType::Iterator PointDataIteratorType;

// Filter
OTB Bot's avatar
STYLE  
OTB Bot committed
50
typedef otb::ImageFileReader<ImageType>       ReaderType;
51
typedef otb::ImageFileWriter<OutputImageType> WriterType;
OTB Bot's avatar
STYLE  
OTB Bot committed
52
typedef otb::ImageFileWriter<ImageType>       WriterInputType;
53

54
OutputImageType::Pointer sift(ImageType::Pointer input,
55 56 57 58 59
                              const unsigned int octaves,
                              const unsigned int scales,
                              const float threshold,
                              const float ratio,
                              const char* siftFileName)
60
{
OTB Bot's avatar
STYLE  
OTB Bot committed
61 62
  typedef otb::ImageToSIFTKeyPointSetFilter<ImageType, PointSetType> SiftFilterType;
  typedef itk::PointSetToImageFilter<PointSetType, OutputImageType>  PointSetFilterType;
63

OTB Bot's avatar
STYLE  
OTB Bot committed
64
  SiftFilterType::Pointer     sift = SiftFilterType::New();
65
  PointSetFilterType::Pointer pointSetFilter = PointSetFilterType::New();
66

67
  sift->SetInput(input);
68 69 70 71
  sift->SetOctavesNumber(octaves);
  sift->SetScalesNumber(scales);
  sift->SetDoGThreshold(threshold);
  sift->SetEdgeThreshold(ratio);
72

73 74
  pointSetFilter->SetInput(sift->GetOutput());
  pointSetFilter->SetOutsideValue(0);
75
  pointSetFilter->SetInsideValue(255);
76
  pointSetFilter->SetSize(input->GetLargestPossibleRegion().GetSize());
77
  pointSetFilter->SetSpacing(input->GetSignedSpacing());
78
  pointSetFilter->SetOrigin(input->GetOrigin());
79
  pointSetFilter->Update();
80

81
  WriterType::Pointer writer = WriterType::New();
82
  writer->SetFileName(siftFileName);
83
  writer->SetInput(pointSetFilter->GetOutput());
84
  //writer->Update();
85

86
  return pointSetFilter->GetOutput();
87 88
}

OTB Bot's avatar
STYLE  
OTB Bot committed
89 90
OutputImageType::Pointer ddm(OutputImageType::Pointer input,
                             const char* ddmFileName)
91
{
92
  typedef itk::DanielssonDistanceMapImageFilter <OutputImageType, OutputImageType> DDMFilterType;
93
  DDMFilterType::Pointer ddmFilter = DDMFilterType::New();
94

95 96 97 98 99 100 101
  ddmFilter->SetInput(input);
  ddmFilter->InputIsBinaryOn();
  ddmFilter->Update();

  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName(ddmFileName);
  writer->SetInput(ddmFilter->GetOutput());
102
  //writer->Update();
103

104 105 106
  return ddmFilter->GetOutput();
}

OTB Bot's avatar
STYLE  
OTB Bot committed
107 108
ImageType::Pointer rotate(ImageType::Pointer input,
                          const unsigned int rotation)
109 110 111
{
  typedef itk::AffineTransform<double, 2> TransformType;
  typedef itk::ResampleImageFilter<ImageType, ImageType>
112
  ResampleFilterType;
113
  ResampleFilterType::Pointer resampler = ResampleFilterType::New();
114

OTB Bot's avatar
STYLE  
OTB Bot committed
115
  TransformType::Pointer          transform = TransformType::New();
116 117
  TransformType::OutputVectorType translation1;
  TransformType::OutputVectorType translation2;
118

119
  const ImageType::SpacingType& spacing = input->GetSignedSpacing();
OTB Bot's avatar
STYLE  
OTB Bot committed
120 121
  const ImageType::PointType&   origin  = input->GetOrigin();
  ImageType::SizeType           size = input->GetLargestPossibleRegion().GetSize();
122

123 124
  const double imageCenterX = origin[0] + spacing[0] * size[0] / 2.0;
  const double imageCenterY = origin[1] + spacing[1] * size[1] / 2.0;
125 126 127
  const double degreesToRadians = atan(1.0) / 45.0;
  const double angle = rotation * degreesToRadians;

128 129 130 131
  translation1[0] = -imageCenterX;
  translation1[1] = -imageCenterY;
  translation2[0] = imageCenterX;
  translation2[1] = imageCenterY;
132

OTB Bot's avatar
STYLE  
OTB Bot committed
133 134 135
  transform->Translate(translation1);
  transform->Rotate2D(-angle, false);
  transform->Translate(translation2);
136

OTB Bot's avatar
STYLE  
OTB Bot committed
137 138 139
  resampler->SetOutputOrigin(origin);
  resampler->SetOutputSpacing(spacing);
  resampler->SetSize(size);
140 141 142 143 144
  resampler->SetTransform(transform);
  resampler->SetInput(input);
  resampler->Update();
  return resampler->GetOutput();
}
145

146
OutputImageType::Pointer invRotate(OutputImageType::Pointer input,
147
                                   const unsigned int rotation)
148 149 150
{
  typedef itk::AffineTransform<double, 2> TransformType;
  typedef itk::ResampleImageFilter<OutputImageType, OutputImageType>
151
  ResampleFilterType;
152
  ResampleFilterType::Pointer resampler = ResampleFilterType::New();
153

OTB Bot's avatar
STYLE  
OTB Bot committed
154
  TransformType::Pointer          transform = TransformType::New();
155 156
  TransformType::OutputVectorType translation1;
  TransformType::OutputVectorType translation2;
157

158
  const ImageType::SpacingType& spacing = input->GetSignedSpacing();
OTB Bot's avatar
STYLE  
OTB Bot committed
159 160
  const ImageType::PointType&   origin  = input->GetOrigin();
  ImageType::SizeType           size = input->GetLargestPossibleRegion().GetSize();
161

162 163
  const double imageCenterX = origin[0] + spacing[0] * size[0] / 2.0;
  const double imageCenterY = origin[1] + spacing[1] * size[1] / 2.0;
164

165 166
  const double degreesToRadians = atan(1.0) / 45.0;
  const double angle = rotation * degreesToRadians;
167

168 169 170 171
  translation1[0] = -imageCenterX;
  translation1[1] = -imageCenterY;
  translation2[0] = imageCenterX;
  translation2[1] = imageCenterY;
172

OTB Bot's avatar
STYLE  
OTB Bot committed
173 174 175
  transform->Translate(translation1);
  transform->Rotate2D(angle, false);
  transform->Translate(translation2);
176

OTB Bot's avatar
STYLE  
OTB Bot committed
177 178 179
  resampler->SetOutputOrigin(origin);
  resampler->SetOutputSpacing(spacing);
  resampler->SetSize(size);
180 181 182 183
  resampler->SetTransform(transform);
  resampler->SetInput(input);
  resampler->Update();
  return resampler->GetOutput();
184 185
}

OTB Bot's avatar
STYLE  
OTB Bot committed
186 187
ImageType::Pointer zoom(ImageType::Pointer input,
                        const unsigned int zoomFactor)
188
{
189
  typedef itk::ShrinkImageFilter<ImageType, ImageType> ShrinkFilterType;
190
  ShrinkFilterType::Pointer shrink = ShrinkFilterType::New();
191

192 193 194
  shrink->SetInput(input);
  shrink->SetShrinkFactors(zoomFactor);
  shrink->Update();
195

196 197 198
  return shrink->GetOutput();
}

OTB Bot's avatar
STYLE  
OTB Bot committed
199 200
OutputImageType::Pointer invZoom(OutputImageType::Pointer input,
                                 const unsigned int zoomFactor)
201 202 203
{
  typedef itk::ExpandImageFilter<OutputImageType, OutputImageType> ExpandFilterType;
  ExpandFilterType::Pointer expand = ExpandFilterType::New();
204

205
  expand->SetInput(input);
206 207
  expand->SetExpandFactors(zoomFactor);
  expand->Update();
208

209
  return expand->GetOutput();
210 211
}

212
ImageType::Pointer contrast(ImageType::Pointer input,
213 214
                            const unsigned int contrastMin,
                            const unsigned int contrastMax)
215
{
216
  typedef itk::RescaleIntensityImageFilter<ImageType, ImageType> RescaleFilterType;
OTB Bot's avatar
STYLE  
OTB Bot committed
217
  RescaleFilterType::Pointer rescaler = RescaleFilterType::New();
218

219 220 221
  rescaler->SetInput(input);
  rescaler->SetOutputMinimum(static_cast<RescaleFilterType::OutputPixelType>(contrastMin));
  rescaler->SetOutputMaximum(static_cast<RescaleFilterType::OutputPixelType>(contrastMax));
222 223
  rescaler->Update();
  return rescaler->GetOutput();
224 225
}

OTB Bot's avatar
STYLE  
OTB Bot committed
226
OutputImageType::Pointer invContrast(OutputImageType::Pointer input,
227 228
                                     const unsigned int itkNotUsed(contrastMin),
                                     const unsigned int itkNotUsed(contrastMax))
229
{
230
  return input;
231 232 233
}

void subtract(OutputImageType::Pointer image1,
234 235
              OutputImageType::Pointer image2,
              const char* subtractFileName)
236
{
237
  typedef itk::SubtractImageFilter<OutputImageType, OutputImageType, ImageType> SubtractFilterType;
OTB Bot's avatar
STYLE  
OTB Bot committed
238 239
  typedef itk::MinimumMaximumImageCalculator<ImageType>                         MaximumCalculatorType;
  typedef itk::RescaleIntensityImageFilter<ImageType, OutputImageType>          OutputRescaleFilterType;
240

OTB Bot's avatar
STYLE  
OTB Bot committed
241
  SubtractFilterType::Pointer    subtract = SubtractFilterType::New();
242
  MaximumCalculatorType::Pointer maximumCalculator = MaximumCalculatorType::New();
243

244 245 246
  subtract->SetInput1(image1);
  subtract->SetInput2(image2);
  subtract->Update();
247

248 249 250 251
  OutputRescaleFilterType::Pointer rescaler = OutputRescaleFilterType::New();
  rescaler->SetInput(subtract->GetOutput());
  rescaler->SetOutputMinimum(0);
  rescaler->SetOutputMaximum(255);
252

253
  WriterType::Pointer writer = WriterType::New();
254
  writer->SetFileName(subtractFileName);
255
  writer->SetInput(rescaler->GetOutput());
256
  //writer->Update();
257

258
  maximumCalculator->SetImage(subtract->GetOutput());
259 260 261
  maximumCalculator->Compute();

  std::cout << "Mix(sub)= " << maximumCalculator->GetMinimum() << " ";
262 263 264 265 266
  std::cout << "Max(sub)= " << maximumCalculator->GetMaximum() << std::endl;
}

int otbImageToSIFTKeyPointSetFilterDistanceMap(int argc, char * argv[])
{
267

268
  if (argc < 10)
OTB Bot's avatar
STYLE  
OTB Bot committed
269
    {
270 271
    std::cout << "Missing arguments " << std::endl;
    return EXIT_FAILURE;
OTB Bot's avatar
STYLE  
OTB Bot committed
272
    }
273

274 275
  // Input Image file name
  const char * infname = argv[1];
276
  const char * outfname = argv[10];
277

278 279
  const unsigned int octaves = atoi(argv[2]);
  const unsigned int scales = atoi(argv[3]);
OTB Bot's avatar
STYLE  
OTB Bot committed
280 281
  const float        threshold = atof(argv[4]);
  const float        ratio = atof(argv[5]);
282

OTB Bot's avatar
STYLE  
OTB Bot committed
283
  // Rotation angle [0, 360[
284
  const unsigned int rotation = atoi(argv[6]);
285

286 287
  // Zoom factor
  const unsigned int zoomFactor = atoi(argv[7]);
288

289 290 291
  // contrast factor
  const unsigned int contrastMin = atoi(argv[8]);
  const unsigned int contrastMax = atoi(argv[9]);
292

293
  //redirect cout to a file
294 295
  std::ofstream   file(outfname);
  std::streambuf* strm_buffer = std::cout.rdbuf(); //save the cout to put it
296 297 298
                                             //back later
  std::cout.rdbuf(file.rdbuf());

299 300 301
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName(infname);
  reader->Update();
302

303 304 305 306 307 308 309 310 311
  ImageType::Pointer _rotated =
    rotate(reader->GetOutput(), rotation);
  ImageType::Pointer _zoomed =
    zoom(reader->GetOutput(), zoomFactor);
  ImageType::Pointer _contrasted =
    contrast(reader->GetOutput(), contrastMin, contrastMax);

  ImageType::Pointer _combined1 = zoom(_rotated, zoomFactor);
  ImageType::Pointer _combined2 = contrast(_combined1, contrastMin, contrastMax);
312

313 314
  OutputImageType::Pointer sift_base =
    sift(reader->GetOutput(), octaves, scales, threshold, ratio, "sift_base.png");
315 316

  OutputImageType::Pointer _sift_rotated =
317
    sift(_rotated, octaves, scales, threshold, ratio, "sift_rotated.png");
318

319 320 321 322 323
  OutputImageType::Pointer _sift_zoomed =
    sift(_zoomed, octaves, scales, threshold, ratio, "sift_zoomed.png");

  OutputImageType::Pointer _sift_contrasted =
    sift(_contrasted, octaves, scales, threshold, ratio, "sift_contrasted.png");
324

325 326
  OutputImageType::Pointer _sift_combined =
    sift(_combined2, octaves, scales, threshold, ratio, "sift_combined.png");
327

328 329
  OutputImageType::Pointer sift_rotated =
    invRotate(_sift_rotated, rotation);
330

331 332
  OutputImageType::Pointer sift_zoomed =
    invZoom(_sift_zoomed, zoomFactor);
333

334 335
  OutputImageType::Pointer sift_contrasted =
    invContrast(_sift_contrasted, contrastMin, contrastMax);
336

337 338 339 340
  OutputImageType::Pointer _sift_combined1 =
    invContrast(_sift_combined, contrastMin, contrastMax);
  OutputImageType::Pointer _sift_combined2 =
    invRotate(_sift_combined1, rotation);
341
  OutputImageType::Pointer sift_combined =
342
    invZoom(_sift_combined2, zoomFactor);
343

344 345 346 347 348
  OutputImageType::Pointer ddm_base = ddm(sift_base, "ddm_base.png");
  OutputImageType::Pointer ddm_rotated = ddm(sift_rotated, "ddm_rotated.png");
  OutputImageType::Pointer ddm_contrasted = ddm(sift_contrasted, "ddm_contrasted.png");
  OutputImageType::Pointer ddm_zoomed = ddm(sift_zoomed, "ddm_zoomed.png");
  OutputImageType::Pointer ddm_combined = ddm(sift_combined, "ddm_combined.png");
349

350
  subtract(ddm_base, ddm_rotated,
351
           "subtract_rotated.png");
352

353
  subtract(ddm_base, ddm_zoomed,
354
           "subtract_zoomed.png");
355

356
  subtract(ddm_base, ddm_contrasted,
357
           "subtract_contrasted.png");
358

359
  subtract(ddm_base, ddm_combined,
360
           "subtract_combined.png");
361

362 363 364
  std::cout.rdbuf(strm_buffer);
  file.close();

365 366
  return EXIT_SUCCESS;
}