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