OTB  6.7.0
Orfeo Toolbox
otbPersistentSamplingFilterBase.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbPersistentSamplingFilterBase_hxx
22 #define otbPersistentSamplingFilterBase_hxx
23 
28 #include "otbMacro.h"
29 #include "otbStopwatch.h"
30 #include "itkProgressReporter.h"
31 
32 namespace otb
33 {
34 
35 template<class TInputImage, class TMaskImage>
38  : m_FieldName(std::string("class"))
39  , m_FieldIndex(0)
40  , m_LayerIndex(0)
41  , m_OutLayerName(std::string("output"))
42  , m_OGRLayerCreationOptions()
43  , m_AdditionalFields()
44  , m_InMemoryInputs()
45  , m_InMemoryOutputs()
46 {
47  this->SetNthOutput(0,TInputImage::New());
48 }
49 
50 
51 template<class TInputImage, class TMaskImage>
52 void
55 {
56  this->SetNthInput(1, const_cast<otb::ogr::DataSource *>( vector ));
57 }
58 
59 template<class TInputImage, class TMaskImage>
63 {
64  if (this->GetNumberOfInputs()<2)
65  {
66  return 0;
67  }
68  return static_cast<const otb::ogr::DataSource *>(this->itk::ProcessObject::GetInput(1));
69 }
70 
71 template<class TInputImage, class TMaskImage>
72 void
74 ::SetMask(const TMaskImage* mask)
75 {
76  this->SetNthInput(2, const_cast<TMaskImage *>( mask ));
77 }
78 
79 template<class TInputImage, class TMaskImage>
80 const TMaskImage*
83 {
84  if (this->GetNumberOfInputs()<3)
85  {
86  return 0;
87  }
88  return static_cast<const TMaskImage *>(this->itk::ProcessObject::GetInput(2));
89 }
90 
91 template<class TInputImage, class TMaskImage>
92 void
94 ::SetOGRLayerCreationOptions(const std::vector<std::string> & options)
95 {
96  m_OGRLayerCreationOptions.clear();
97  m_OGRLayerCreationOptions = options;
98 }
99 
100 template<class TInputImage, class TMaskImage>
101 const std::vector<std::string>&
104 {
105  return m_OGRLayerCreationOptions;
106 }
107 
108 template<class TInputImage, class TMaskImage>
109 void
112 {
113  Superclass::GenerateOutputInformation();
114 
115  // Get OGR field index
116  const otb::ogr::DataSource* vectors = this->GetOGRData();
117  otb::ogr::Layer::const_iterator featIt = vectors->GetLayer(m_LayerIndex).begin();
118  int fieldIndex = featIt->ogr().GetFieldIndex(this->m_FieldName.c_str());
119  if (fieldIndex < 0)
120  {
121  itkGenericExceptionMacro("Field named "<<this->m_FieldName<<" not found!");
122  }
123  this->m_FieldIndex = fieldIndex;
124 
125  const MaskImageType *mask = this->GetMask();
126  if (mask)
127  {
128  const InputImageType *input = this->GetInput();
129  if (mask->GetLargestPossibleRegion() !=
130  input->GetLargestPossibleRegion() )
131  {
132  itkGenericExceptionMacro("Mask and input image have a different size!");
133  }
134  if (mask->GetOrigin() != input->GetOrigin())
135  {
136  itkGenericExceptionMacro("Mask and input image have a different origin!");
137  }
138  if (mask->GetSignedSpacing() != input->GetSignedSpacing())
139  {
140  itkGenericExceptionMacro("Mask and input image have a different spacing!");
141  }
142  }
143 }
144 
145 template<class TInputImage, class TMaskImage>
146 void
149 {
150  InputImageType *input = const_cast<InputImageType*>(this->GetInput());
151  MaskImageType *mask = const_cast<MaskImageType*>(this->GetMask());
152 
153  RegionType requested = this->GetOutput()->GetRequestedRegion();
154  RegionType emptyRegion = input->GetLargestPossibleRegion();
155  emptyRegion.SetSize(0,0);
156  emptyRegion.SetSize(1,0);
157 
158  input->SetRequestedRegion(emptyRegion);
159 
160  if (mask)
161  {
162  mask->SetRequestedRegion(requested);
163  }
164 }
165 
166 template <class TInputImage, class TMaskImage>
167 void
170 {
171  this->AllocateOutputs();
172  this->BeforeThreadedGenerateData();
173 
174  // Split the data into in-memory layers
175  this->DispatchInputVectors();
176 
177  // struct to store filter pointer
178  VectorThreadStruct str;
179  str.Filter = this;
180 
181  // Get the output pointer
182  //const InputImageType *outputPtr = this->GetOutput();
183 
184  this->GetMultiThreader()->SetNumberOfThreads( this->GetNumberOfThreads() );
185  this->GetMultiThreader()->SetSingleMethod(this->VectorThreaderCallback, &str);
186 
187  // multithread the execution
188  this->GetMultiThreader()->SingleMethodExecute();
189 
190  // gather the data from in-memory output layers
191  this->GatherOutputVectors();
192 
193  this->AfterThreadedGenerateData();
194 }
195 
196 template <class TInputImage, class TMaskImage>
197 void
200 {
201  Superclass::AllocateOutputs();
202 
203  ogr::DataSource* vectors = const_cast<ogr::DataSource*>(this->GetOGRData());
204  ogr::Layer inLayer = vectors->GetLayer(m_LayerIndex);
205 
206  unsigned int numberOfThreads = this->GetNumberOfThreads();
207 
208  // Prepare temporary input
209  this->m_InMemoryInputs.clear();
210  this->m_InMemoryInputs.reserve(numberOfThreads);
211  std::string tmpLayerName("thread");
212  OGRSpatialReference * oSRS = nullptr;
213  if (inLayer.GetSpatialRef())
214  {
215  oSRS = inLayer.GetSpatialRef()->Clone();
216  }
217  OGRFeatureDefn &layerDefn = inLayer.GetLayerDefn();
218  //std::vector<ogr::Layer> tmpLayers;
219  for (unsigned int i=0 ; i < numberOfThreads ; i++)
220  {
222  ogr::Layer tmpLayer = tmpOgrDS->CreateLayer(
223  tmpLayerName,
224  oSRS,
225  inLayer.GetGeomType());
226  // add field definitions
227  for (int k=0 ; k < layerDefn.GetFieldCount() ; k++)
228  {
229  OGRFieldDefn originDefn(layerDefn.GetFieldDefn(k));
230  ogr::FieldDefn fieldDefn(originDefn);
231  tmpLayer.CreateField(fieldDefn);
232  }
233  this->m_InMemoryInputs.push_back(tmpOgrDS);
234  }
235 
236  // Prepare in-memory outputs
237  this->m_InMemoryOutputs.clear();
238  this->m_InMemoryOutputs.reserve(numberOfThreads);
239  tmpLayerName = std::string("threadOut");
240  for (unsigned int i=0 ; i < numberOfThreads ; i++)
241  {
242  std::vector<OGRDataPointer> tmpContainer;
243  // iterate over outputs, only process ogr::DataSource
244  for (unsigned int k=0 ; k < this->GetNumberOfOutputs() ; k++)
245  {
246  ogr::DataSource* realOutput = dynamic_cast<ogr::DataSource *>(
248  if (realOutput)
249  {
250  ogr::Layer realLayer = realOutput->GetLayersCount() == 1
251  ? realOutput->GetLayer(0)
252  : realOutput->GetLayer(m_OutLayerName);
253  OGRFeatureDefn &outLayerDefn = realLayer.GetLayerDefn();
255  ogr::Layer tmpLayer = tmpOutput->CreateLayer(
256  tmpLayerName, oSRS, realLayer.GetGeomType());
257  // add field definitions
258  for (int f=0 ; f < outLayerDefn.GetFieldCount() ; f++)
259  {
260  OGRFieldDefn originDefn(outLayerDefn.GetFieldDefn(f));
261  tmpLayer.CreateField(originDefn);
262  }
263  tmpContainer.push_back(tmpOutput);
264  }
265  }
266  this->m_InMemoryOutputs.push_back(tmpContainer);
267  }
268 
269  if (oSRS)
270  {
271  oSRS->Release();
272  }
273 }
274 
275 template <class TInputImage, class TMaskImage>
276 void
279 {
280  // clean temporary inputs
281  this->m_InMemoryInputs.clear();
282 
283  // gather temporary outputs and write to output
284  const otb::ogr::DataSource* vectors = this->GetOGRData();
286  unsigned int count = 0;
287  for (unsigned int k=0 ; k < this->GetNumberOfOutputs() ; k++)
288  {
289  ogr::DataSource* realOutput = dynamic_cast<ogr::DataSource *>(
291  if (realOutput)
292  {
293  this->FillOneOutput(count, realOutput, bool(vectors == realOutput));
294  count++;
295  }
296  }
297 
298  chrono.Stop();
299  otbMsgDebugMacro(<< "Writing OGR points took " << chrono.GetElapsedMilliseconds() << " ms");
300  this->m_InMemoryOutputs.clear();
301 }
302 
303 template <class TInputImage, class TMaskImage>
304 void
306 ::FillOneOutput(unsigned int outIdx, ogr::DataSource* outDS, bool update)
307 {
308  ogr::Layer outLayer = outDS->GetLayersCount() == 1
309  ? outDS->GetLayer(0)
310  : outDS->GetLayer(m_OutLayerName);
311 
312  OGRErr err = outLayer.ogr().StartTransaction();
313  if (err != OGRERR_NONE)
314  {
315  itkExceptionMacro(<< "Unable to start transaction for OGR layer " << outLayer.ogr().GetName() << ".");
316  }
317 
318  unsigned int numberOfThreads = this->GetNumberOfThreads();
319  for (unsigned int thread=0 ; thread < numberOfThreads ; thread++)
320  {
321  ogr::Layer inLayer = this->m_InMemoryOutputs[thread][outIdx]->GetLayerChecked(0);
322  if (!inLayer)
323  {
324  continue;
325  }
326 
327  ogr::Layer::const_iterator tmpIt = inLayer.begin();
328  // This test only uses 1 input, not compatible with multiple OGRData inputs
329  if (update)
330  {
331  // Update mode
332  for(; tmpIt!=inLayer.end(); ++tmpIt)
333  {
334  outLayer.SetFeature( *tmpIt );
335  }
336  }
337  else
338  {
339  // Copy mode
340  for(; tmpIt!=inLayer.end(); ++tmpIt)
341  {
342  ogr::Feature dstFeature(outLayer.GetLayerDefn());
343  dstFeature.SetFrom( *tmpIt, TRUE );
344  outLayer.CreateFeature( dstFeature );
345  }
346  }
347  }
348 
349  err = outLayer.ogr().CommitTransaction();
350  if (err != OGRERR_NONE)
351  {
352  itkExceptionMacro(<< "Unable to commit transaction for OGR layer " << outLayer.ogr().GetName() << ".");
353  }
354 }
355 
356 template <class TInputImage, class TMaskImage>
357 void
360 {
361  // Retrieve inputs
362  TInputImage* inputImage = const_cast<TInputImage*>(this->GetInput());
363  TInputImage* outputImage = this->GetOutput();
364  RegionType requestedRegion = outputImage->GetRequestedRegion();
365 
366  itk::ProgressReporter progress( this, threadid, layerForThread.GetFeatureCount(true) );
367 
368  // Loop across the features in the layer (filtered by requested region in BeforeTGD already)
369  ogr::Layer::const_iterator featIt = layerForThread.begin();
370  for(; featIt!=layerForThread.end(); ++featIt)
371  {
372  // Compute the intersection of thread region and polygon bounding region, called "considered region"
373  // This need not be done in ThreadedGenerateData and could be pre-processed and cached before filter execution if needed
374  RegionType consideredRegion = FeatureBoundingRegion(inputImage, featIt);
375  bool regionNotEmpty = consideredRegion.Crop(requestedRegion);
376  if (regionNotEmpty)
377  {
378  this->PrepareFeature(*featIt,threadid);
379  this->ExploreGeometry(*featIt, featIt->ogr().GetGeometryRef(),consideredRegion,threadid);
380  }
381  progress.CompletedPixel();
382  }
383 }
384 
385 template <class TInputImage, class TMaskImage>
386 void
389  OGRGeometry* geom,
390  RegionType& region,
391  itk::ThreadIdType& threadid)
392 {
393  typename TInputImage::PointType imgPoint;
394  typename TInputImage::IndexType imgIndex;
395 
396  switch (geom->getGeometryType())
397  {
398  case wkbPoint:
399  case wkbPoint25D:
400  {
401  OGRPoint* castPoint = dynamic_cast<OGRPoint*>(geom);
402  if (castPoint == nullptr) break;
403 
404  imgPoint[0] = castPoint->getX();
405  imgPoint[1] = castPoint->getY();
406  const TInputImage* img = this->GetInput();
407  const TMaskImage* mask = this->GetMask();
408  img->TransformPhysicalPointToIndex(imgPoint,imgIndex);
409  if ((mask == nullptr) || mask->GetPixel(imgIndex))
410  {
411  this->ProcessSample(feature, imgIndex, imgPoint, threadid);
412  }
413  break;
414  }
415  case wkbLineString:
416  case wkbLineString25D:
417  {
418  OGRLineString* castLineString = dynamic_cast<OGRLineString*>(geom);
419 
420  if (castLineString == nullptr) break;
421  this->ProcessLine(feature,castLineString,region,threadid);
422  break;
423  }
424  case wkbPolygon:
425  case wkbPolygon25D:
426  {
427  OGRPolygon* castPolygon = dynamic_cast<OGRPolygon*>(geom);
428  if (castPolygon == nullptr) break;
429  this->ProcessPolygon(feature,castPolygon,region,threadid);
430  break;
431  }
432  case wkbMultiPoint:
433  case wkbMultiPoint25D:
434  case wkbMultiLineString:
435  case wkbMultiLineString25D:
436  case wkbMultiPolygon:
437  case wkbMultiPolygon25D:
438  case wkbGeometryCollection:
439  case wkbGeometryCollection25D:
440  {
441  OGRGeometryCollection *geomCollection = dynamic_cast<OGRGeometryCollection*>(geom);
442  if (geomCollection)
443  {
444  unsigned int nbGeom = geomCollection->getNumGeometries();
445  for (unsigned int i=0 ; i < nbGeom ; ++i)
446  {
447  this->ExploreGeometry(feature,
448  geomCollection->getGeometryRef(i),
449  region,
450  threadid);
451  }
452  }
453  else
454  {
455  otbWarningMacro("Geometry not recognized as a collection : " << geom->getGeometryName());
456  }
457  break;
458  }
459  default:
460  {
461  otbWarningMacro("Geometry not handled: " << geom->getGeometryName());
462  break;
463  }
464  }
465 }
466 
467 template <class TInputImage, class TMaskImage>
468 void
471  OGRLineString* line,
472  RegionType& region,
473  itk::ThreadIdType& threadid)
474 {
475  OGRPolygon tmpPolygon;
476  OGRLinearRing ring;
477  ring.addPoint(0.0,0.0,0.0);
478  ring.addPoint(1.0,0.0,0.0);
479  ring.addPoint(1.0,1.0,0.0);
480  ring.addPoint(0.0,1.0,0.0);
481  ring.addPoint(0.0,0.0,0.0);
482  tmpPolygon.addRing(&ring);
483  const TInputImage* img = this->GetInput();
484  TMaskImage* mask = const_cast<TMaskImage*>(this->GetMask());
485  typename TInputImage::IndexType imgIndex;
486  typename TInputImage::PointType imgPoint;
487  typename TInputImage::SpacingType imgAbsSpacing = img->GetSignedSpacing();
488  if (imgAbsSpacing[0] < 0) imgAbsSpacing[0] = -imgAbsSpacing[0];
489  if (imgAbsSpacing[1] < 0) imgAbsSpacing[1] = -imgAbsSpacing[1];
490 
491  if (mask)
492  {
493  // For pixels in consideredRegion and not masked
494  typedef MaskedIteratorDecorator<
496  itk::ImageRegionConstIterator<TMaskImage> > MaskedIteratorType;
497  MaskedIteratorType it(mask, mask, region);
498  it.GoToBegin();
499  while (!it.IsAtEnd())
500  {
501  imgIndex = it.GetIndex();
502  img->TransformIndexToPhysicalPoint(imgIndex,imgPoint);
503  bool isInside = this->IsSampleOnLine(line,imgPoint,imgAbsSpacing,tmpPolygon);
504  if (isInside)
505  {
506  this->ProcessSample(feature,imgIndex, imgPoint, threadid);
507  }
508  ++it;
509  }
510  }
511  else
512  {
514  NoValueIteratorType it(img,region);
515  it.GoToBegin();
516  while (!it.IsAtEnd())
517  {
518  imgIndex = it.GetIndex();
519  img->TransformIndexToPhysicalPoint(imgIndex,imgPoint);
520  bool isInside = this->IsSampleOnLine(line,imgPoint,imgAbsSpacing,tmpPolygon);
521  if (isInside)
522  {
523  this->ProcessSample(feature,imgIndex, imgPoint, threadid);
524  }
525  ++it;
526  }
527  }
528 }
529 
530 template <class TInputImage, class TMaskImage>
531 void
534  OGRPolygon* polygon,
535  RegionType& region,
536  itk::ThreadIdType& threadid)
537 {
538  const TInputImage* img = this->GetInput();
539  TMaskImage* mask = const_cast<TMaskImage*>(this->GetMask());
540  typename TInputImage::IndexType imgIndex;
541  typename TInputImage::PointType imgPoint;
542  OGRPoint tmpPoint;
543 
544  if (mask)
545  {
546  // For pixels in consideredRegion and not masked
547  typedef MaskedIteratorDecorator<
549  itk::ImageRegionConstIterator<TMaskImage> > MaskedIteratorType;
550  MaskedIteratorType it(mask, mask, region);
551  it.GoToBegin();
552  while (!it.IsAtEnd())
553  {
554  imgIndex = it.GetIndex();
555  img->TransformIndexToPhysicalPoint(imgIndex,imgPoint);
556  tmpPoint.setX(imgPoint[0]);
557  tmpPoint.setY(imgPoint[1]);
558  bool isInside = this->IsSampleInsidePolygon(polygon,&tmpPoint);
559  if (isInside)
560  {
561  this->ProcessSample(feature,imgIndex, imgPoint, threadid);
562  }
563  ++it;
564  }
565  }
566  else
567  {
569  NoValueIteratorType it(img,region);
570  it.GoToBegin();
571  while (!it.IsAtEnd())
572  {
573  imgIndex = it.GetIndex();
574  img->TransformIndexToPhysicalPoint(imgIndex,imgPoint);
575  tmpPoint.setX(imgPoint[0]);
576  tmpPoint.setY(imgPoint[1]);
577  bool isInside = this->IsSampleInsidePolygon(polygon,&tmpPoint);
578  if (isInside)
579  {
580  this->ProcessSample(feature,imgIndex, imgPoint, threadid);
581  }
582  ++it;
583  }
584  }
585 }
586 
587 template <class TInputImage, class TMaskImage>
588 void
591  typename TInputImage::IndexType& ,
592  typename TInputImage::PointType& ,
594 {
595  itkExceptionMacro("Method ProcessSample not implemented !");
596 }
597 
598 template <class TInputImage, class TMaskImage>
599 void
603 {
604  // Nothing to do here
605 }
606 
607 template <class TInputImage, class TMaskImage>
608 inline bool
610 ::IsSampleInsidePolygon(OGRPolygon* poly,
611  OGRPoint* tmpPoint)
612 {
613  bool ret = poly->getExteriorRing()->isPointInRing(tmpPoint);
614  if (ret)
615  {
616  for (int k=0 ; k<poly->getNumInteriorRings() ; k++)
617  {
618  if (poly->getInteriorRing(k)->isPointInRing(tmpPoint))
619  {
620  ret = false;
621  break;
622  }
623  }
624  }
625  return ret;
626 }
627 
628 template <class TInputImage, class TMaskImage>
629 inline bool
631 ::IsSampleOnLine(OGRLineString* line,
632  typename TInputImage::PointType& position,
633  typename TInputImage::SpacingType& absSpacing,
634  OGRPolygon& tmpPolygon)
635 {
636  tmpPolygon.getExteriorRing()->setPoint(0
637  ,position[0]-0.5*absSpacing[0]
638  ,position[1]-0.5*absSpacing[1]
639  ,0.0);
640  tmpPolygon.getExteriorRing()->setPoint(1
641  ,position[0]+0.5*absSpacing[0]
642  ,position[1]-0.5*absSpacing[1]
643  ,0.0);
644  tmpPolygon.getExteriorRing()->setPoint(2
645  ,position[0]+0.5*absSpacing[0]
646  ,position[1]+0.5*absSpacing[1]
647  ,0.0);
648  tmpPolygon.getExteriorRing()->setPoint(3
649  ,position[0]-0.5*absSpacing[0]
650  ,position[1]+0.5*absSpacing[1]
651  ,0.0);
652  tmpPolygon.getExteriorRing()->setPoint(4
653  ,position[0]-0.5*absSpacing[0]
654  ,position[1]-0.5*absSpacing[1]
655  ,0.0);
656  return line->Intersects(&tmpPolygon);
657 }
658 
659 template<class TInputImage, class TMaskImage>
662 ::FeatureBoundingRegion(const TInputImage* image, otb::ogr::Layer::const_iterator& featIt) const
663 {
664  // otb::ogr wrapper is incomplete and leaky abstraction is inevitable here
665  OGREnvelope envelope;
666  featIt->GetGeometry()->getEnvelope(&envelope);
667  itk::Point<double, 2> lowerPoint, upperPoint;
668  lowerPoint[0] = envelope.MinX;
669  lowerPoint[1] = envelope.MinY;
670  upperPoint[0] = envelope.MaxX;
671  upperPoint[1] = envelope.MaxY;
672 
673  typename TInputImage::IndexType lowerIndex;
674  typename TInputImage::IndexType upperIndex;
675 
676  image->TransformPhysicalPointToIndex(lowerPoint, lowerIndex);
677  image->TransformPhysicalPointToIndex(upperPoint, upperIndex);
678 
679  // swap coordinate to keep lowerIndex as start index
680  if (lowerIndex[0] > upperIndex[0])
681  {
682  int tmp = lowerIndex[0];
683  lowerIndex[0] = upperIndex[0];
684  upperIndex[0] = tmp;
685  }
686  if (lowerIndex[1] > upperIndex[1])
687  {
688  int tmp = lowerIndex[1];
689  lowerIndex[1] = upperIndex[1];
690  upperIndex[1] = tmp;
691  }
692 
693  RegionType region;
694  region.SetIndex(lowerIndex);
695  region.SetUpperIndex(upperIndex);
696 
697  return region;
698 }
699 
700 template<class TInputImage, class TMaskImage>
701 void
704 {
705  TInputImage* outputImage = this->GetOutput();
706  ogr::DataSource* vectors = const_cast<ogr::DataSource*>(this->GetOGRData());
707  ogr::Layer inLayer = vectors->GetLayer(m_LayerIndex);
708 
709  const RegionType& requestedRegion = outputImage->GetRequestedRegion();
710  itk::ContinuousIndex<double> startIndex(requestedRegion.GetIndex());
711  itk::ContinuousIndex<double> endIndex(requestedRegion.GetUpperIndex());
712  startIndex[0] += -0.5;
713  startIndex[1] += -0.5;
714  endIndex[0] += 0.5;
715  endIndex[1] += 0.5;
716  itk::Point<double, 2> startPoint;
717  itk::Point<double, 2> endPoint;
718  outputImage->TransformContinuousIndexToPhysicalPoint(startIndex, startPoint);
719  outputImage->TransformContinuousIndexToPhysicalPoint(endIndex, endPoint);
720 
721  // create geometric extent
722  OGRPolygon tmpPolygon;
723  OGRLinearRing ring;
724  ring.addPoint(startPoint[0],startPoint[1],0.0);
725  ring.addPoint(startPoint[0],endPoint[1] ,0.0);
726  ring.addPoint(endPoint[0] ,endPoint[1] ,0.0);
727  ring.addPoint(endPoint[0] ,startPoint[1],0.0);
728  ring.addPoint(startPoint[0],startPoint[1],0.0);
729  tmpPolygon.addRing(&ring);
730 
731  inLayer.SetSpatialFilter(&tmpPolygon);
732 
733  unsigned int numberOfThreads = this->GetNumberOfThreads();
734  std::vector<ogr::Layer> tmpLayers;
735  tmpLayers.reserve(numberOfThreads);
736  for (unsigned int i=0 ; i<numberOfThreads ; i++)
737  {
738  tmpLayers.push_back(this->GetInMemoryInput(i));
739  }
740 
741  const unsigned int nbFeatThread = std::ceil(inLayer.GetFeatureCount(true) / (float) numberOfThreads);
742  //assert(nbFeatThread > 0);
743 
744  OGRFeatureDefn &layerDefn = inLayer.GetLayerDefn();
745  ogr::Layer::const_iterator featIt = inLayer.begin();
746  unsigned int counter=0;
747  unsigned int cptFeat = 0;
748  for(; featIt!=inLayer.end(); ++featIt)
749  {
750  ogr::Feature dstFeature(layerDefn);
751  dstFeature.SetFrom( *featIt, TRUE );
752  dstFeature.SetFID(featIt->GetFID());
753  tmpLayers[counter].CreateFeature( dstFeature );
754  cptFeat++;
755  if (cptFeat > nbFeatThread && (counter + 1) < numberOfThreads)
756  {
757  counter++;
758  cptFeat=0;
759  }
760  }
761 
762  inLayer.SetSpatialFilter(nullptr);
763 }
764 
765 template<class TInputImage, class TMaskImage>
766 void
769 {
770  TInputImage *inputImage = const_cast<TInputImage*>(this->GetInput());
771  inputImage->UpdateOutputInformation();
772 
773  ogr::Layer inLayer = inputDS->GetLayer(this->GetLayerIndex());
774 
775  bool updateMode = false;
776  if (inputDS == outputDS)
777  {
778  updateMode = true;
779  // Check m_OutLayerName is same as input layer name
780  m_OutLayerName = inLayer.GetName();
781  }
782 
783  // First get list of current fields
784  OGRFeatureDefn &layerDefn = inLayer.GetLayerDefn();
785  std::map<std::string, OGRFieldType> currentFields;
786  for (int k=0 ; k<layerDefn.GetFieldCount() ; k++)
787  {
788  OGRFieldDefn fieldDefn(layerDefn.GetFieldDefn(k));
789  std::string currentName(fieldDefn.GetNameRef());
790  currentFields[currentName] = fieldDefn.GetType();
791  }
792 
793  ogr::Layer outLayer = inLayer;
794  if (!updateMode)
795  {
796  std::string projectionRefWkt = this->GetInput()->GetProjectionRef();
797  bool projectionInformationAvailable = !projectionRefWkt.empty();
798  OGRSpatialReference * oSRS = nullptr;
799  if(projectionInformationAvailable)
800  {
801  oSRS = static_cast<OGRSpatialReference *>(OSRNewSpatialReference(projectionRefWkt.c_str()));
802  }
803  // Create layer
804  outLayer = outputDS->CreateLayer(
805  this->GetOutLayerName(),
806  oSRS,
807  wkbPoint,
808  this->GetOGRLayerCreationOptions());
809  // Copy existing fields
810  for (int k=0 ; k<layerDefn.GetFieldCount() ; k++)
811  {
812  OGRFieldDefn fieldDefn(layerDefn.GetFieldDefn(k));
813  outLayer.CreateField(fieldDefn);
814  }
815 
816  if (oSRS)
817  {
818  oSRS->Release();
819  }
820  }
821 
822  // Add new fields
823  for (unsigned int k=0 ; k<m_AdditionalFields.size() ; k++)
824  {
825  OGRFieldDefn ogrFieldDefinition(m_AdditionalFields[k].Name.c_str(),m_AdditionalFields[k].Type);
826  ogrFieldDefinition.SetWidth( m_AdditionalFields[k].Width );
827  ogrFieldDefinition.SetPrecision( m_AdditionalFields[k].Precision );
828  ogr::FieldDefn fieldDef(ogrFieldDefinition);
829  // test if field is already present
830  if (currentFields.count(fieldDef.GetName()))
831  {
832  // test the field type
833  if (currentFields[fieldDef.GetName()] != fieldDef.GetType())
834  {
835  itkExceptionMacro("Field name "<< fieldDef.GetName() << " already exists with a different type!");
836  }
837  }
838  else
839  {
840  outLayer.CreateField(fieldDef);
841  }
842  }
843 }
844 
845 
846 template<class TInputImage, class TMaskImage>
847 void
850 {
851  this->m_AdditionalFields.clear();
852 }
853 
854 template<class TInputImage, class TMaskImage>
855 void
857 ::CreateAdditionalField(std::string name,
858  OGRFieldType type,
859  int width,
860  int precision)
861 {
862  SimpleFieldDefn defn;
863  defn.Name = name;
864  defn.Type = type;
865  defn.Width = width;
866  defn.Precision = precision;
867  this->m_AdditionalFields.push_back(defn);
868 }
869 
870 template<class TInputImage, class TMaskImage>
871 const std::vector<
873  ::SimpleFieldDefn>&
876 {
877  return this->m_AdditionalFields;
878 }
879 
880 template<class TInputImage, class TMaskImage>
884 {
885  VectorThreadStruct *str = (VectorThreadStruct*)(((itk::MultiThreader::ThreadInfoStruct *)(arg))->UserData);
886 
887  int threadId = ((itk::MultiThreader::ThreadInfoStruct *)(arg))->ThreadID;
888  int threadCount = ((itk::MultiThreader::ThreadInfoStruct *)(arg))->NumberOfThreads;
889 
890  ogr::Layer layer = str->Filter->GetInMemoryInput(threadId);
891 
892  if (threadId < threadCount)
893  {
894  str->Filter->ThreadedGenerateVectorData(layer,threadId);
895  }
896 
898 }
899 
900 template<class TInputImage, class TMaskImage>
903 ::GetInMemoryInput(unsigned int threadId)
904 {
905  if (threadId >= m_InMemoryInputs.size())
906  {
907  itkExceptionMacro(<< "Requested in-memory input layer not available " << threadId << " (total size : "<< m_InMemoryInputs.size() <<").");
908  }
909  return m_InMemoryInputs[threadId]->GetLayerChecked(0);
910 }
911 
912 template<class TInputImage, class TMaskImage>
915 ::GetInMemoryOutput(unsigned int threadId, unsigned int index)
916 {
917  if (threadId >= m_InMemoryOutputs.size())
918  {
919  itkExceptionMacro(<< "Requested in-memory output layer not available " << threadId << " (total size : "<< m_InMemoryOutputs.size() <<").");
920  }
921  if (index >= m_InMemoryOutputs[threadId].size())
922  {
923  itkExceptionMacro(<< "Requested output dataset not available " << index << " (available : "<< m_InMemoryOutputs[threadId].size() <<").");
924  }
925  return m_InMemoryOutputs[threadId][index]->GetLayerChecked(0);
926 }
927 
928 } // end namespace otb
929 
930 #endif
virtual void ProcessPolygon(const ogr::Feature &feature, OGRPolygon *polygon, RegionType &region, itk::ThreadIdType &threadid)
virtual void PrepareFeature(const ogr::Feature &feature, itk::ThreadIdType &threadid)
Base class for persistent filter doing sampling tasks.
Collection of geometric objects.
DurationType GetElapsedMilliseconds() const
#define ITK_THREAD_RETURN_VALUE
virtual void InitializeOutputDataSource(ogr::DataSource *inputDS, ogr::DataSource *outputDS)
void CreateField(FieldDefn const &field, bool bApproxOK=true)
OGRFeatureDefn & GetLayerDefn() const
static Stopwatch StartNew()
Layer GetLayer(vcl_size_t i)
ogr::Layer GetInMemoryOutput(unsigned int threadId, unsigned int index=0)
Layer of geometric objets.
bool IsSampleOnLine(OGRLineString *line, typename TInputImage::PointType &position, typename TInputImage::SpacingType &absSpacing, OGRPolygon &tmpPolygon)
Stopwatch timer.
Definition: otbStopwatch.h:41
RegionType FeatureBoundingRegion(const TInputImage *image, otb::ogr::Layer::const_iterator &featIt) const
void SetOGRData(const ogr::DataSource *vector)
virtual void ThreadedGenerateVectorData(const ogr::Layer &layerForThread, itk::ThreadIdType threadid)
#define ITK_THREAD_RETURN_TYPE
Implementation class for Feature iterator. This iterator is a single pass iterator. We may fetch the Feature referenced by an iterator previously stored, but never resume the iteration after a call to Layer::begin(), Layer::start_at(), Layer::CreateFeature(), Layer::DeleteFeature(), Layer::GetFeature(), Layer::SetFeature(), nor fork the iteration.
virtual void FillOneOutput(unsigned int outIdx, ogr::DataSource *outDS, bool update)
const_iterator begin() const
itk::Index< Monteverdi_DIMENSION > IndexType
Definition: mvdTypes.h:133
void CreateFeature(Feature feature)
#define otbMsgDebugMacro(x)
Definition: otbMacro.h:64
Decorate an iterator to ignore masked pixels.
OGRwkbGeometryType GetGeomType() const
static Pointer New()
OGRSpatialReference const * GetSpatialRef() const
int GetFeatureCount(bool doForceComputation) const
void CreateAdditionalField(std::string name, OGRFieldType type, int width=0, int precision=0)
#define otbWarningMacro(x)
Definition: otbMacro.h:67
std::string GetName() const
Field name accessor.
virtual void ProcessLine(const ogr::Feature &feature, OGRLineString *line, RegionType &region, itk::ThreadIdType &threadid)
const std::vector< SimpleFieldDefn > & GetAdditionalFields()
unsigned int ThreadIdType
DataObject * GetInput(const DataObjectIdentifierType &key)
std::string GetProjectionRef() const
OGRLayer & ogr()
void SetOGRLayerCreationOptions(const std::vector< std::string > &options)
int GetLayersCount() const
void SetFrom(Feature const &rhs, int *map, bool mustForgive=true)
const std::vector< std::string > & GetOGRLayerCreationOptions()
Geometric objet with descriptive fields.
VectorImageType::SpacingType SpacingType
Definition: mvdTypes.h:181
virtual void ProcessSample(const ogr::Feature &feature, typename TInputImage::IndexType &imgIndex, typename TInputImage::PointType &imgPoint, itk::ThreadIdType &threadid)
static ITK_THREAD_RETURN_TYPE VectorThreaderCallback(void *arg)
bool IsSampleInsidePolygon(OGRPolygon *poly, OGRPoint *tmpPoint)
Layer CreateLayer(std::string const &name, OGRSpatialReference *poSpatialRef=nullptr, OGRwkbGeometryType eGType=wkbUnknown, std::vector< std::string > const &papszOptions=std::vector< std::string >())
void SetFeature(Feature feature)
VectorImageType::PointType PointType
Definition: mvdTypes.h:189
std::string GetName() const
Encapsulation of OGRFieldDefn: field definition.
void ExploreGeometry(const ogr::Feature &feature, OGRGeometry *geom, RegionType &region, itk::ThreadIdType &threadid)
ogr::Layer GetInMemoryInput(unsigned int threadId)
void SetSpatialFilter(OGRGeometry const *spatialFilter)
const_iterator end() const
OGRFieldType GetType() const
Field type accessor.
DataObject * GetOutput(const DataObjectIdentifierType &key)