otbHomologousPointsExtraction.cxx 21.5 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 25

#include "otbWrapperApplicationFactory.h"

// Elevation handler
#include "otbWrapperElevationParametersHandler.h"

26
#ifdef OTB_USE_SIFTFAST
27
#include "otbSiftFastImageFilter.h"
28
#else
29
#include "otbImageToSIFTKeyPointSetFilter.h"
30
#endif
31 32 33 34 35
#include "otbImageToSURFKeyPointSetFilter.h"
#include "otbKeyPointSetsMatchingFilter.h"
#include "otbMultiToMonoChannelExtractROI.h"
#include "otbKeyPointSetsMatchingFilter.h"

36
#include "otbGenericRSTransform.h"
37 38
#include "otbOGRDataSourceWrapper.h"
#include "ogrsf_frmts.h"
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56

namespace otb
{
namespace Wrapper
{

class HomologousPointsExtraction: public Application
{
public:
  /** Standard class typedefs. */
  typedef HomologousPointsExtraction      Self;
  typedef Application                     Superclass;
  typedef itk::SmartPointer<Self>         Pointer;
  typedef itk::SmartPointer<const Self>   ConstPointer;

  typedef FloatVectorImageType::PixelType                    RealVectorType;
  typedef FloatVectorImageType::InternalPixelType            ValueType;
  typedef itk::PointSet<RealVectorType,2>                    PointSetType;
57 58 59

  // This is ugly ...
  #ifdef OTB_USE_SIFTFAST
60
  typedef SiftFastImageFilter<FloatImageType,PointSetType>   SiftFilterType;
61
  #else
Julien Michel's avatar
Julien Michel committed
62
  typedef ImageToSIFTKeyPointSetFilter<FloatImageType,PointSetType> SiftFilterType;
63 64
  #endif

65 66
  typedef ImageToSURFKeyPointSetFilter<FloatImageType,PointSetType> SurfFilterType;

67
  typedef itk::Statistics::EuclideanDistanceMetric<RealVectorType> DistanceType;
68 69 70 71 72 73 74
  typedef otb::KeyPointSetsMatchingFilter<PointSetType,
                                          DistanceType>      MatchingFilterType;

  typedef MatchingFilterType::LandmarkListType               LandmarkListType;
  typedef PointSetType::PointType                            PointType;
  typedef otb::MultiToMonoChannelExtractROI<ValueType,
                                            ValueType>      ExtractChannelFilterType;
75
  typedef otb::GenericRSTransform<>                         RSTransformType;
76 77 78 79 80 81 82

  /** Standard macro */
  itkNewMacro(Self);

  itkTypeMacro(HomologousPointsExtraction, otb::Wrapper::Application);

private:
83
  void DoInit() override
84 85
  {
    SetName("HomologousPointsExtraction");
Julien Michel's avatar
Julien Michel committed
86
    SetDocName("Homologous points extraction");
87 88
    SetDescription("Compute homologous points between images using keypoints");
    SetDocLongDescription("This application allows computing homologous points between images using keypoints. "
89
      " SIFT or SURF keypoints can be used and the band on which keypoints are computed can be set independently for both images."
90 91 92
      " The application offers two modes :"
      " the first is the full mode where keypoints are extracted from the full extent of both images"
      " (please note that in this mode large image file are not supported). "
93
      "The second mode, called geobins, allows one to set-up spatial binning to get fewer points"
94
      " spread across the entire image. "
95 96 97 98
      "In this mode, the corresponding spatial bin in the second image is estimated using geographical"
      " transform or sensor modelling, and is padded according to the user defined precision. Last, in"
      " both modes the application can filter matches whose colocalisation in first image exceed this precision. "
      "The elevation parameters are to deal more precisely with sensor modelling in case of sensor geometry data. "
99
      "The outvector option allows creating a vector file with segments corresponding to the localisation error between the matches."
100 101 102
      " It can be useful to assess the precision of a registration for instance."
      " The vector file is always reprojected to EPSG:4326 to allow display in a GIS."
      " This is done via reprojection or by applying the image sensor models.");
103 104
    // Documentation
    SetDocName("Homologous Points Extraction");
105
    SetDocLimitations("Full mode does not handle large images.");
Julien Michel's avatar
Julien Michel committed
106
    SetDocSeeAlso("RefineSensorModel");
107 108
    SetDocAuthors("OTB-Team");

109
    AddDocTag(Tags::FeatureExtraction);
110 111 112 113

    AddParameter(ParameterType_InputImage, "in1", "Input Image 1");
    SetParameterDescription("in1"," First input image");

114 115 116 117 118
    AddParameter(ParameterType_Int,"band1","Input band 1");
    SetParameterDescription("band1","Index of the band from input image 1 to use for keypoints extraction");
    SetMinimumParameterIntValue("band1",1);
    SetDefaultParameterInt("band1",1);

119 120 121
    AddParameter(ParameterType_InputImage, "in2", "Input Image 2");
    SetParameterDescription("in2"," Second input image");

122 123 124 125 126
    AddParameter(ParameterType_Int,"band2","Input band 2");
    SetParameterDescription("band2","Index of the band from input image 1 to use for keypoints extraction");
    SetMinimumParameterIntValue("band2",1);
    SetDefaultParameterInt("band2",1);

127 128 129 130 131 132 133 134
    AddParameter(ParameterType_OutputFilename, "out", "Output file with tie points");
    SetParameterDescription("out", "File containing the list of tie points");

    AddParameter(ParameterType_OutputFilename, "outvector", "Output vector file with tie points");
    SetParameterDescription("outvector", "File containing segments representing matches ");
    MandatoryOff("outvector");
    DisableParameter("outvector");

135 136
    AddParameter(ParameterType_Choice,"algorithm","Keypoints detection algorithm");
    SetParameterDescription("algorithm","Choice of the detection algorithm to use");
137

138
    AddChoice("algorithm.surf","SURF algorithm");
139
    AddChoice("algorithm.sift","SIFT algorithm");
140

141 142 143 144 145
    AddParameter(ParameterType_Float,"threshold","Distance threshold for matching");
    SetParameterDescription("threshold","The distance threshold for matching.");
    SetMinimumParameterFloatValue("threshold",0.0);
    SetDefaultParameterFloat("threshold",0.6);

146
    AddParameter(ParameterType_Bool,"backmatching","Use back-matching to filter matches.");
147 148
    SetParameterDescription("backmatching","If set to true, matches should be consistent in both ways.");

149
    AddParameter(ParameterType_Choice,"mode","Keypoints search mode");
150

151 152 153
    AddChoice("mode.full","Extract and match all keypoints (no streaming)");
    SetParameterDescription("mode.full","Extract and match all keypoints, loading both images entirely into memory");

154
    AddChoice("mode.geobins","Search keypoints in small spatial bins regularly spread across first image");
155
    SetParameterDescription("mode.geobins","This method allows retrieving a set of tie points regulary spread across image 1. Corresponding bins in image 2 are retrieved using sensor and geographical information if available. The first bin position takes into account the margin parameter. Bins are cropped to the largest image region shrunk by the margin parameter for both in1 and in2 images.");
156

157 158 159 160 161
    AddParameter(ParameterType_Int,"mode.geobins.binsize","Size of bin");
    SetParameterDescription("mode.geobins.binsize","Radius of the spatial bin in pixels");
    SetDefaultParameterInt("mode.geobins.binsize",256);
    SetMinimumParameterIntValue("mode.geobins.binsize",1);

162 163 164 165
    AddParameter(ParameterType_Int,"mode.geobins.binsizey","Size of bin (y direction)");
    SetParameterDescription("mode.geobins.binsizey","Radius of the spatial bin in pixels (y direction). If not set, the mode.geobins.binsize value is used.");
    SetMinimumParameterIntValue("mode.geobins.binsizey",1);
    MandatoryOff("mode.geobins.binsizey");
166

167 168 169 170 171
    AddParameter(ParameterType_Int,"mode.geobins.binstep","Steps between bins");
    SetParameterDescription("mode.geobins.binstep","Steps between bins in pixels");
    SetDefaultParameterInt("mode.geobins.binstep",256);
    SetMinimumParameterIntValue("mode.geobins.binstep",1);

172 173 174 175 176
    AddParameter(ParameterType_Int,"mode.geobins.binstepy","Steps between bins (y direction)");
    SetParameterDescription("mode.geobins.binstepy","Steps between bins in pixels (y direction). If not set, the mode.geobins.binstep value is used.");
    SetMinimumParameterIntValue("mode.geobins.binstepy",1);
    MandatoryOff("mode.geobins.binstepy");

177 178 179 180 181
    AddParameter(ParameterType_Int,"mode.geobins.margin","Margin from image border to start/end bins (in pixels)");
    SetParameterDescription("mode.geobins.margin","Margin from image border to start/end bins (in pixels)");
    SetMinimumParameterIntValue("mode.geobins.margin",0);
    SetDefaultParameterInt("mode.geobins.margin",10);

182 183 184
    AddParameter(ParameterType_Float,"precision","Estimated precision of the colocalisation function (in pixels).");
    SetParameterDescription("precision","Estimated precision of the colocalisation function in pixels");
    SetDefaultParameterFloat("precision",0.);
185

186
    AddParameter(ParameterType_Bool,"mfilter","Filter points according to geographical or sensor based colocalisation");
187
    SetParameterDescription("mfilter","If enabled, this option allows one to filter matches according to colocalisation from sensor or geographical information, using the given tolerancy expressed in pixels");
188

189
    AddParameter(ParameterType_Bool,"2wgs84","If enabled, points from second image will be exported in WGS84");
190

191 192 193
    // Elevation
    ElevationParametersHandler::AddElevationParameters(this, "elev");

Julien Michel's avatar
Julien Michel committed
194
    // Doc example parameter settings
195 196
    SetDocExampleParameterValue("in1", "sensor_stereo_left.tif");
    SetDocExampleParameterValue("in2","sensor_stereo_right.tif");
Julien Michel's avatar
Julien Michel committed
197 198
    SetDocExampleParameterValue("mode","full");
    SetDocExampleParameterValue("out","homologous.txt");
199

200
    SetOfficialDocLink();
201 202
  }

203
  void DoUpdateParameters() override
204 205 206 207
  {

  }

208
  void Match(FloatImageType * im1, FloatImageType * im2, RSTransformType * rsTransform, RSTransformType * rsTransform1ToWGS84,RSTransformType * rsTransform2ToWGS84, std::ofstream & file, OGRMultiLineString * mls = nullptr)
209
  {
210
    MatchingFilterType::Pointer matchingFilter = MatchingFilterType::New();
211

212 213
    if(GetParameterString("algorithm")=="sift")
      {
214
      otbAppLogINFO("Using SIFT points");
215 216
      SiftFilterType::Pointer sift1 = SiftFilterType::New();
      sift1->SetInput(im1);
217

218 219
      SiftFilterType::Pointer sift2 = SiftFilterType::New();
      sift2->SetInput(im2);
220

221
      sift1->Update();
222

223
      otbAppLogINFO("Found " << sift1->GetOutput()->GetNumberOfPoints()<<" sift points in image 1.");
224

225
      sift2->Update();
226

227
      otbAppLogINFO("Found " << sift2->GetOutput()->GetNumberOfPoints()<<" sift points in image 2.");
228

229 230 231 232 233 234 235
      matchingFilter->SetInput1(sift1->GetOutput());
      matchingFilter->SetInput2(sift2->GetOutput());
      }
    else if(GetParameterString("algorithm")=="surf")
      {
      SurfFilterType::Pointer surf1 = SurfFilterType::New();
      surf1->SetInput(im1);
236

237 238
      SurfFilterType::Pointer surf2 = SurfFilterType::New();
      surf2->SetInput(im2);
239

240
      otbAppLogINFO("Doing update");
241
      surf1->Update();
242

243
      otbAppLogINFO("Found " << surf1->GetOutput()->GetNumberOfPoints()<<" surf points in image 1.");
244

245
      surf2->Update();
246

247
      otbAppLogINFO("Found " << surf2->GetOutput()->GetNumberOfPoints()<<" surf points in image 2.");
248

249 250
      matchingFilter->SetInput1(surf1->GetOutput());
      matchingFilter->SetInput2(surf2->GetOutput());
251
      matchingFilter->SetDistanceThreshold(GetParameterFloat("threshold"));
252
      matchingFilter->SetUseBackMatching(GetParameterInt("backmatching"));
253
      }
254

255
    try
256 257
      {

258
      matchingFilter->Update();
259

260
      LandmarkListType::Pointer landmarks = matchingFilter->GetOutput();
261

262
      otbAppLogINFO("Found " << landmarks->Size() <<" homologous points.");
263

264
      unsigned int discarded  = 0;
265

266 267
      for (LandmarkListType::Iterator it = landmarks->Begin();
           it != landmarks->End(); ++it)
268
        {
269 270
        PointType point1 = it.Get()->GetPoint1();
        PointType point2 = it.Get()->GetPoint2();
271

272
        double error = 0;
273
        PointType pprime1,pprime2;
274

275
        bool filtered = false;
276

277
        if(GetParameterInt("mfilter"))
278
          {
279
          pprime1 = rsTransform->TransformPoint(point1);
280
          error = std::sqrt((point2[0]-pprime1[0])*(point2[0]-pprime1[0])+(point2[1]-pprime1[1])*(point2[1]-pprime1[1]));
281

282
          if(error>GetParameterFloat("precision")*std::sqrt(std::abs(im2->GetSignedSpacing()[0]*im2->GetSignedSpacing()[1])))
283 284 285
            {
            filtered = true;
            }
286
          }
287

288
        if(!filtered)
289
          {
290
          if(GetParameterInt("2wgs84"))
291
            {
292
            pprime2 = rsTransform2ToWGS84->TransformPoint(point2);
293

294
            file<<point1[0]<<"\t"<<point1[1]<<"\t"<<pprime2[0]<<"\t"<<pprime2[1]<<std::endl;
295 296 297 298 299
            }
          else
            {
            file<<point1[0]<<"\t"<<point1[1]<<"\t"<<point2[0]<<"\t"<<point2[1]<<std::endl;
            }
300

301 302 303 304 305 306 307 308 309 310 311
          if(mls)
            {
             pprime1 = rsTransform1ToWGS84->TransformPoint(point1);
             pprime2 = rsTransform2ToWGS84->TransformPoint(point2);

             OGRLineString ls;
             ls.addPoint(pprime1[0],pprime1[1]);
             ls.addPoint(pprime2[0],pprime2[1]);
             mls->addGeometry(&ls);
            }

312 313 314
          }
        else
          {
315
          ++discarded;
316 317
          }
        }
318 319
      otbAppLogINFO(<<discarded<<" points discarded");
      }
320
    catch(itk::ExceptionObject &)
321 322
      {
      // silent catch
323 324
      }
  }
325

326

327
  void DoExecute() override
328
  {
329 330
    OGRMultiLineString mls;

331 332 333
    FloatVectorImageType::Pointer image1 = this->GetParameterImage("in1");
    FloatVectorImageType::Pointer image2 = this->GetParameterImage("in2");

334
    // Setting up RS Transform
335
    RSTransformType::Pointer rsTransform = RSTransformType::New();
336 337 338 339
    rsTransform->SetInputKeywordList(image1->GetImageKeywordlist());
    rsTransform->SetInputProjectionRef(image1->GetProjectionRef());
    rsTransform->SetOutputKeywordList(image2->GetImageKeywordlist());
    rsTransform->SetOutputProjectionRef(image2->GetProjectionRef());
340

341
    RSTransformType::Pointer rsTransform1ToWGS84 = RSTransformType::New();
342 343
    rsTransform1ToWGS84->SetInputKeywordList(image1->GetImageKeywordlist());
    rsTransform1ToWGS84->SetInputProjectionRef(image1->GetProjectionRef());
344

345
    RSTransformType::Pointer rsTransform2ToWGS84 = RSTransformType::New();
346 347
    rsTransform2ToWGS84->SetInputKeywordList(image2->GetImageKeywordlist());
    rsTransform2ToWGS84->SetInputProjectionRef(image2->GetProjectionRef());
348

349 350
    // Setting up output file
    std::ofstream file;
351
    file.open(GetParameterString("out"));
352 353 354 355 356
    file<<std::fixed;
    file.precision(12);

    // Setting up channel extractors
    ExtractChannelFilterType::Pointer extractChannel1 = ExtractChannelFilterType::New();
357
    extractChannel1->SetInput(image1);
358
    extractChannel1->SetChannel(GetParameterInt("band1"));
359

360
    ExtractChannelFilterType::Pointer extractChannel2 = ExtractChannelFilterType::New();
361
    extractChannel2->SetInput(image2);
362
    extractChannel2->SetChannel(GetParameterInt("band2"));
363

364 365
    // Setup the DEM Handler
    otb::Wrapper::ElevationParametersHandler::SetupDEMHandlerFromElevationParameters(this,"elev");
366

367 368 369
    rsTransform->InstantiateTransform();
    rsTransform1ToWGS84->InstantiateTransform();
    rsTransform2ToWGS84->InstantiateTransform();
370

371

372 373 374
    if(GetParameterString("mode")=="full")
      {
      // Launch detection on whole images
375
      Match(extractChannel1->GetOutput(),extractChannel2->GetOutput(),rsTransform,rsTransform1ToWGS84,rsTransform2ToWGS84,file,&mls);
376

377 378
      }
    else if(GetParameterString("mode")=="geobins")
379
      {
380
      // Compute binning on first image
381
      FloatImageType::SizeType size = image1->GetLargestPossibleRegion().GetSize();
382 383 384
      unsigned int bin_size_x = GetParameterInt("mode.geobins.binsize");
      unsigned int bin_size_y = bin_size_x;

385 386
      unsigned int image_border_margin = GetParameterInt("mode.geobins.margin");

387 388 389 390
      if(IsParameterEnabled("mode.geobins.binsizey"))
        {
        bin_size_y = GetParameterInt("mode.geobins.binsizey");
        }
391

392 393 394 395 396 397 398
      unsigned int bin_step_x = GetParameterInt("mode.geobins.binstep");
      unsigned int bin_step_y = bin_step_x;

      if(IsParameterEnabled("mode.geobins.binstepy"))
        {
        bin_step_y = GetParameterInt("mode.geobins.binstepy");
        }
399

400 401
      unsigned int nb_bins_x = static_cast<unsigned int>(std::ceil(static_cast<float>(size[0]-2*image_border_margin)/(bin_size_x + bin_step_x)));
      unsigned int nb_bins_y = static_cast<unsigned int>(std::ceil(static_cast<float>(size[1]-2*image_border_margin)/(bin_size_y + bin_step_y)));
402

OTB Bot's avatar
STYLE  
OTB Bot committed
403
      for(unsigned int i = 0; i<nb_bins_x; ++i)
404
        {
405 406
        for(unsigned int j = 0; j<nb_bins_y; ++j)
          {
407 408
          unsigned int startx = image_border_margin + i*(bin_size_x + bin_step_x);
          unsigned int starty = image_border_margin + j*(bin_size_y + bin_step_y);
409

410

411 412 413
          FloatImageType::SizeType size1;
          FloatImageType::IndexType index1;
          FloatImageType::RegionType region1;
414

415 416
          index1[0]=startx;
          index1[1]=starty;
417 418
          size1[0] = bin_size_x;
          size1[1] = bin_size_y;
419

420 421
          region1.SetIndex(index1);
          region1.SetSize(size1);
422

423
          FloatImageType::RegionType largestRegion = image1->GetLargestPossibleRegion();
424 425 426

          largestRegion.ShrinkByRadius(image_border_margin);
          region1.Crop(largestRegion);
427

428
          otbAppLogINFO("("<<i+1<<"/"<<nb_bins_x<<", "<<j+1<<"/"<<nb_bins_y<<") Considering region1 : "<<region1.GetIndex()<<", "<<region1.GetSize());
429

430

431
          extractChannel1->SetExtractionRegion(region1);
432 433


434
          // We need to find the corresponding region in image 2
435
          itk::ContinuousIndex<double,2> ul_i, ur_i, lr_i, ll_i;
436 437
          FloatImageType::PointType ul1, ur1, ll1, lr1, p1, p2, p3, p4;
          itk::ContinuousIndex<double,2> i1, i2, i3, i4, i_min, i_max;
438

439 440 441 442 443 444 445 446 447 448 449 450 451 452
          ul_i[0] = static_cast<double>(startx) - 0.5;
          ul_i[1] = static_cast<double>(starty) - 0.5;
          ur_i = ul_i;
          lr_i = ul_i;
          ll_i = ul_i;
          ur_i[0] += static_cast<double>(bin_size_x);
          lr_i[0] += static_cast<double>(bin_size_x);
          lr_i[1] += static_cast<double>(bin_size_y);
          ll_i[1] += static_cast<double>(bin_size_y);

          image1->TransformContinuousIndexToPhysicalPoint(ul_i,ul1);
          image1->TransformContinuousIndexToPhysicalPoint(ur_i,ur1);
          image1->TransformContinuousIndexToPhysicalPoint(lr_i,lr1);
          image1->TransformContinuousIndexToPhysicalPoint(ll_i,ll1);
453

454 455 456 457
          p1 = rsTransform->TransformPoint(ul1);
          p2 = rsTransform->TransformPoint(ur1);
          p3 = rsTransform->TransformPoint(lr1);
          p4 = rsTransform->TransformPoint(ll1);
458

459 460 461 462
          image2->TransformPhysicalPointToContinuousIndex(p1,i1);
          image2->TransformPhysicalPointToContinuousIndex(p2,i2);
          image2->TransformPhysicalPointToContinuousIndex(p3,i3);
          image2->TransformPhysicalPointToContinuousIndex(p4,i4);
463

464 465
          i_min[0] = std::min(std::min(i1[0],i2[0]),std::min(i3[0],i4[0]));
          i_min[1] = std::min(std::min(i1[1],i2[1]),std::min(i3[1],i4[1]));
466

467 468
          i_max[0] = std::max(std::max(i1[0],i2[0]),std::max(i3[0],i4[0]));
          i_max[1] = std::max(std::max(i1[1],i2[1]),std::max(i3[1],i4[1]));
469

470 471
          FloatImageType::IndexType index2;
          FloatImageType::SizeType size2;
472

473 474
          index2[0] = std::floor(i_min[0]);
          index2[1] = std::floor(i_min[1]);
475

476 477
          size2[0] = std::ceil(i_max[0]-i_min[0]);
          size2[1] = std::ceil(i_max[1]-i_min[1]);
478

479 480 481 482
          FloatImageType::RegionType region2;
          region2.SetIndex(index2);
          region2.SetSize(size2);
          region2.PadByRadius(static_cast<unsigned int>(GetParameterInt("precision")));
483

484
          FloatImageType::RegionType largestRegion2 = image2->GetLargestPossibleRegion();
485 486 487
          largestRegion2.ShrinkByRadius(image_border_margin);

          if(region2.Crop(largestRegion2))
OTB Bot's avatar
STYLE  
OTB Bot committed
488
            {
489
            otbAppLogINFO("Corresponding region2 is "<<region2.GetIndex()<<", "<<region2.GetSize());
490

491
            extractChannel2->SetExtractionRegion(region2);
492

493
            Match(extractChannel1->GetOutput(),extractChannel2->GetOutput(),rsTransform,rsTransform1ToWGS84,rsTransform2ToWGS84,file,&mls);
494
            }
495 496 497 498
          else
            {
            otbAppLogINFO("Corresponding region2 is out of range: "<<region2.GetIndex()<<", "<<region2.GetSize());
            }
499
          }
OTB Bot's avatar
STYLE  
OTB Bot committed
500
        }
501
      }
502
      file.close();
503 504 505 506

      if(IsParameterEnabled("outvector"))
        {
        // Create the datasource (for matches export)
507
        otb::ogr::Layer layer(nullptr, false);
508
        otb::ogr::DataSource::Pointer ogrDS;
509

510 511 512
        ogrDS = otb::ogr::DataSource::New(GetParameterString("outvector"), otb::ogr::DataSource::Modes::Overwrite);
        std::string projref = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
        OGRSpatialReference oSRS(projref.c_str());
513

514 515 516 517
        // and create the layer
        layer = ogrDS->CreateLayer("matches", &oSRS, wkbMultiLineString);
        OGRFeatureDefn & defn = layer.GetLayerDefn();
        ogr::Feature feature(defn);
518

519 520 521
        feature.SetGeometry(&mls);
        layer.CreateFeature(feature);
        }
522 523 524 525 526
  }
};
}
}

527

528
  OTB_APPLICATION_EXPORT(otb::Wrapper::HomologousPointsExtraction)