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