OTB  6.7.0
Orfeo Toolbox
otbListSampleGenerator.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 otbListSampleGenerator_hxx
22 #define otbListSampleGenerator_hxx
23 
24 #include "otbListSampleGenerator.h"
25 
27 #include "otbRemoteSensingRegion.h"
29 
30 #include "otbMacro.h"
31 
32 namespace otb
33 {
34 
35 /*template <class TVectorData>
36 void printVectorData(TVectorData * vectorData, string msg = "")
37 {
38  typedef TVectorData VectorDataType;
39  typedef itk::PreOrderTreeIterator<typename VectorDataType::DataTreeType> TreeIteratorType;
40 
41  TreeIteratorType itVector(vectorData->GetDataTree());
42  itVector.GoToBegin();
43 
44  if (!msg.empty())
45  {
46  std::cout<< msg << std::endl;
47  }
48 
49  while (!itVector.IsAtEnd())
50  {
51  if (itVector.Get()->IsPolygonFeature())
52  {
53  std::cout << itVector.Get()->GetNodeTypeAsString() << std::endl;
54  for (unsigned int itPoints = 0; itPoints < itVector.Get()->GetPolygonExteriorRing()->GetVertexList()->Size(); itPoints++)
55  {
56  std::cout << "vertex[" << itPoints << "]: " << itVector.Get()->GetPolygonExteriorRing()->GetVertexList()->GetElement(itPoints) <<std::endl;
57  }
58  std::cout << "Polygon bounding region:\n" << itVector.Get()->GetPolygonExteriorRing()->GetBoundingRegion() << std::endl;
59  }
60  ++itVector;
61  }
62 }*/
63 
64 template<class TImage, class TVectorData>
67  m_MaxTrainingSize(-1),
68  m_MaxValidationSize(-1),
69  m_ValidationTrainingProportion(0.0),
70  m_BoundByMin(true),
71  m_PolygonEdgeInclusion(false),
72  m_NumberOfClasses(0),
73  m_ClassKey("Class"),
74  m_ClassMinSize(-1)
75 {
76  this->SetNumberOfRequiredInputs(2);
77  this->SetNumberOfRequiredOutputs(4);
78 
79  // Register the outputs
80  this->itk::ProcessObject::SetNthOutput(0, this->MakeOutput(0).GetPointer());
81  this->itk::ProcessObject::SetNthOutput(1, this->MakeOutput(1).GetPointer());
82  this->itk::ProcessObject::SetNthOutput(2, this->MakeOutput(2).GetPointer());
83  this->itk::ProcessObject::SetNthOutput(3, this->MakeOutput(3).GetPointer());
84 
85  m_RandomGenerator = RandomGeneratorType::GetInstance();
86 }
87 
88 template <class TImage, class TVectorData>
89 void
91 ::SetInput(const ImageType * image)
92 {
93  this->ProcessObject::SetNthInput(0, const_cast<ImageType *>(image));
94 }
95 
96 template <class TImage, class TVectorData>
97 const TImage *
99 ::GetInput() const
100 {
101  if (this->GetNumberOfInputs() < 1)
102  {
103  return nullptr;
104  }
105 
106  return static_cast<const ImageType *>(this->ProcessObject::GetInput(0));
107 }
108 
109 template <class TImage, class TVectorData>
110 void
113 {
114  this->ProcessObject::SetNthInput(1, const_cast<VectorDataType *>(vectorData));
115 
116  //printVectorData(vectorData);
117 
118 }
119 
120 template <class TImage, class TVectorData>
121 const TVectorData *
124 {
125  if (this->GetNumberOfInputs() < 2)
126  {
127  return nullptr;
128  }
129 
130  return static_cast<const VectorDataType *>(this->ProcessObject::GetInput(1));
131 }
132 
133 template <class TImage, class TVectorData>
137 {
138  DataObjectPointer output;
139  switch (idx)
140  {
141  case 0:
142  output = static_cast<itk::DataObject*>(ListSampleType::New().GetPointer());
143  break;
144  case 1:
145  output = static_cast<itk::DataObject*>(ListLabelType::New().GetPointer());
146  break;
147  case 2:
148  output = static_cast<itk::DataObject*>(ListSampleType::New().GetPointer());
149  break;
150  case 3:
151  output = static_cast<itk::DataObject*>(ListLabelType::New().GetPointer());
152  break;
153  default:
154  output = static_cast<itk::DataObject*>(ListSampleType::New().GetPointer());
155  break;
156  }
157  return output;
158 }
159 // Get the Training ListSample
160 template <class TImage, class TVectorData>
164 {
165  return dynamic_cast<ListSampleType*>(this->itk::ProcessObject::GetOutput(0));
166 }
167 // Get the Training label ListSample
168 template <class TImage, class TVectorData>
172 {
173  return dynamic_cast<ListLabelType*>(this->itk::ProcessObject::GetOutput(1));
174 }
175 
176 // Get the validation ListSample
177 template <class TImage, class TVectorData>
181 {
182  return dynamic_cast<ListSampleType*>(this->itk::ProcessObject::GetOutput(2));
183 }
184 
185 
186 // Get the validation label ListSample
187 template <class TImage, class TVectorData>
191 {
192  return dynamic_cast<ListLabelType*>(this->itk::ProcessObject::GetOutput(3));
193 }
194 
195 template <class TImage, class TVectorData>
196 void
199 {
200  ImagePointerType img = static_cast<ImageType *>(this->ProcessObject::GetInput(0));
201 
202  if(img.IsNotNull())
203  {
204 
205  // Requested regions will be generated during GenerateData
206  // call. For now request an empty region so as to avoid requesting
207  // the largest possible region (fixes bug #943 )
208  typename ImageType::RegionType dummyRegion;
209  typename ImageType::SizeType dummySize;
210  dummySize.Fill(0);
211  dummyRegion.SetSize(dummySize);
212  img->SetRequestedRegion(dummyRegion);
213  }
214 }
215 
216 
217 template <class TImage, class TVectorData>
218 void
221 {
222  // Get the inputs
223  ImagePointerType image = const_cast<ImageType*>(this->GetInput());
224  VectorDataPointerType vectorData = const_cast<VectorDataType*>(this->GetInputVectorData());
225 
226  // Get the outputs
227  ListSamplePointerType trainingListSample = this->GetTrainingListSample();
228  ListLabelPointerType trainingListLabel = this->GetTrainingListLabel();
229  ListSamplePointerType validationListSample = this->GetValidationListSample();
230  ListLabelPointerType validationListLabel = this->GetValidationListLabel();
231 
232  // Gather some information about the relative size of the classes
233  // We would like to have the same number of samples per class
234  this->GenerateClassStatistics();
235 
236  this->ComputeClassSelectionProbability();
237 
238  // Clear the sample lists
239  trainingListSample->Clear();
240  trainingListLabel->Clear();
241  validationListSample->Clear();
242  validationListLabel->Clear();
243 
244  // Set MeasurementVectorSize for each sample list
245  trainingListSample->SetMeasurementVectorSize(image->GetNumberOfComponentsPerPixel());
246  // stores label as integers,so put the size to 1
247  trainingListLabel->SetMeasurementVectorSize(1);
248  validationListSample->SetMeasurementVectorSize(image->GetNumberOfComponentsPerPixel());
249  // stores label as integers,so put the size to 1
250  validationListLabel->SetMeasurementVectorSize(1);
251 
252  m_ClassesSamplesNumberTraining.clear();
253  m_ClassesSamplesNumberValidation.clear();
254 
255  typename ImageType::RegionType imageLargestRegion = image->GetLargestPossibleRegion();
256 
257  TreeIteratorType itVector(vectorData->GetDataTree());
258  for (itVector.GoToBegin(); !itVector.IsAtEnd(); ++itVector)
259  {
260  if (itVector.Get()->IsPolygonFeature())
261  {
262  PolygonPointerType exteriorRing = itVector.Get()->GetPolygonExteriorRing();
263 
264  typename ImageType::RegionType polygonRegion =
265  otb::TransformPhysicalRegionToIndexRegion(exteriorRing->GetBoundingRegion(),
266  image.GetPointer());
267 
268  const bool hasIntersection = polygonRegion.Crop(imageLargestRegion);
269  if (!hasIntersection)
270  {
271  continue;
272  }
273 
274  image->SetRequestedRegion(polygonRegion);
275  image->PropagateRequestedRegion();
276  image->UpdateOutputData();
277 
279  IteratorType it(image, polygonRegion);
280 
281  for (it.GoToBegin(); !it.IsAtEnd(); ++it)
282  {
284  image->TransformIndexToPhysicalPoint(it.GetIndex(), point);
285 
286  if ( exteriorRing->IsInside(point) ||
287  (this->GetPolygonEdgeInclusion() && exteriorRing->IsOnEdge(point)) )
288  {
289  PolygonListPointerType interiorRings = itVector.Get()->GetPolygonInteriorRings();
290 
291  bool isInsideInteriorRing = false;
292  for (typename PolygonListType::Iterator interiorRing = interiorRings->Begin();
293  interiorRing != interiorRings->End();
294  ++interiorRing)
295  {
296  if ( interiorRing.Get()->IsInside(point)
297  || (this->GetPolygonEdgeInclusion() && interiorRing.Get()->IsOnEdge(point)) )
298  {
299  isInsideInteriorRing = true;
300  break;
301  }
302  }
303  if (isInsideInteriorRing)
304  {
305  continue; // skip this pixel and continue
306  }
307 
308  double randomValue = m_RandomGenerator->GetUniformVariate(0.0, 1.0);
309  if (randomValue < m_ClassesProbTraining[itVector.Get()->GetFieldAsInt(m_ClassKey)])
310  {
311  //Add the sample to the training list
312  trainingListSample->PushBack(it.Get());
313  trainingListLabel->PushBack(itVector.Get()->GetFieldAsInt(m_ClassKey));
314  m_ClassesSamplesNumberTraining[itVector.Get()->GetFieldAsInt(m_ClassKey)] += 1;
315  }
316  else if (randomValue < m_ClassesProbTraining[itVector.Get()->GetFieldAsInt(m_ClassKey)]
317  + m_ClassesProbValidation[itVector.Get()->GetFieldAsInt(m_ClassKey)])
318  {
319  //Add the sample to the validation list
320  validationListSample->PushBack(it.Get());
321  validationListLabel->PushBack(itVector.Get()->GetFieldAsInt(m_ClassKey));
322  m_ClassesSamplesNumberValidation[itVector.Get()->GetFieldAsInt(m_ClassKey)] += 1;
323  }
324  //Note: some samples may not be used at all
325  }
326  }
327  }
328  }
329 
330  assert(trainingListSample->Size() == trainingListLabel->Size());
331  assert(validationListSample->Size() == validationListLabel->Size());
332  this->UpdateProgress(1.0f);
333 }
334 
335 template <class TImage, class TVectorData>
336 void
339 {
340  m_ClassesSize.clear();
341 
342  ImageType* image = const_cast<ImageType*> (this->GetInput());
343  typename VectorDataType::ConstPointer vectorData = this->GetInputVectorData();
344 
345  // Compute cumulative area of all polygons of each class
346  TreeIteratorType itVector(vectorData->GetDataTree());
347  for (itVector.GoToBegin(); !itVector.IsAtEnd(); ++itVector)
348  {
349  DataNodeType* datanode = itVector.Get();
350  if (datanode->IsPolygonFeature())
351  {
352  double area = GetPolygonAreaInPixelsUnits(datanode, image);
353  m_ClassesSize[datanode->GetFieldAsInt(m_ClassKey)] += area;
354  }
355  }
356  m_NumberOfClasses = m_ClassesSize.size();
357 }
358 
359 template <class TImage, class TVectorData>
360 void
363 {
364  m_ClassesProbTraining.clear();
365  m_ClassesProbValidation.clear();
366 
367  // Sanity check
368  if (m_ClassesSize.empty())
369  {
370  itkGenericExceptionMacro(<< "No training sample found inside image");
371  }
372 
373  // Go through the classes size to find the smallest one
374  double minSize = itk::NumericTraits<double>::max();
375  for (std::map<ClassLabelType, double>::const_iterator itmap = m_ClassesSize.begin();
376  itmap != m_ClassesSize.end();
377  ++itmap)
378  {
379  if (minSize > itmap->second)
380  {
381  minSize = itmap->second;
382  }
383  }
384 
385  // Apply the proportion between training and validation samples (all training by default)
386  double minSizeTraining = minSize * (1.0 - m_ValidationTrainingProportion);
387  double minSizeValidation = minSize * m_ValidationTrainingProportion;
388 
389  // Apply the limit if specified by the user
390  if(m_BoundByMin)
391  {
392  if ((m_MaxTrainingSize != -1) && (m_MaxTrainingSize < minSizeTraining))
393  {
394  minSizeTraining = m_MaxTrainingSize;
395  }
396  if ((m_MaxValidationSize != -1) && (m_MaxValidationSize < minSizeValidation))
397  {
398  minSizeValidation = m_MaxValidationSize;
399  }
400  }
401  // Compute the probability selection for each class
402  for (std::map<ClassLabelType, double>::const_iterator itmap = m_ClassesSize.begin();
403  itmap != m_ClassesSize.end();
404  ++itmap)
405  {
406  if (m_BoundByMin)
407  {
408  m_ClassesProbTraining[itmap->first] = minSizeTraining / itmap->second;
409  m_ClassesProbValidation[itmap->first] = minSizeValidation / itmap->second;
410  }
411  else
412  {
413  long int maxSizeT = (itmap->second)*(1.0 - m_ValidationTrainingProportion);
414  long int maxSizeV = (itmap->second)*m_ValidationTrainingProportion;
415 
416  // Check if max sizes respect the maximum bounds
417  double correctionRatioTrain = 1.0;
418  if((m_MaxTrainingSize > -1) && (m_MaxTrainingSize < maxSizeT))
419  {
420  correctionRatioTrain = (double)(m_MaxTrainingSize) / (double)(maxSizeT);
421  }
422  double correctionRatioValid = 1.0;
423  if((m_MaxValidationSize > -1) && (m_MaxValidationSize < maxSizeV))
424  {
425  correctionRatioValid = (double)(m_MaxValidationSize) / (double)(maxSizeV);
426  }
427  double correctionRatio = std::min(correctionRatioTrain,correctionRatioValid);
428  m_ClassesProbTraining[itmap->first] = correctionRatio*(1.0 - m_ValidationTrainingProportion);
429  m_ClassesProbValidation[itmap->first] = correctionRatio*m_ValidationTrainingProportion;
430  }
431  }
432 }
433 template <class TImage, class TVectorData>
434 double
437 {
438  const double pixelArea = std::abs(image->GetSignedSpacing()[0] * image->GetSignedSpacing()[1]);
439 
440  // Compute area of exterior ring in pixels
441  PolygonPointerType exteriorRing = polygonDataNode->GetPolygonExteriorRing();
442  double area = exteriorRing->GetArea() / pixelArea;
443 
444  // Remove contribution of all interior rings
445  PolygonListPointerType interiorRings = polygonDataNode->GetPolygonInteriorRings();
446  for (typename PolygonListType::Iterator interiorRing = interiorRings->Begin();
447  interiorRing != interiorRings->End();
448  ++interiorRing)
449  {
450  area -= interiorRing.Get()->GetArea() / pixelArea;
451  }
452 
453  return area;
454 }
455 
456 template <class TImage, class TVectorData>
457 void
459 ::PrintSelf(std::ostream& os, itk::Indent indent) const
460 {
461  os << indent << "* MaxTrainingSize: " << m_MaxTrainingSize << "\n";
462  os << indent << "* MaxValidationSize: " << m_MaxValidationSize << "\n";
463  os << indent << "* Proportion: " << m_ValidationTrainingProportion << "\n";
464  os << indent << "* Input data:\n";
465  if (m_ClassesSize.empty())
466  {
467  os << indent << "Empty\n";
468  }
469  else
470  {
471  for (std::map<ClassLabelType, double>::const_iterator itmap = m_ClassesSize.begin();
472  itmap != m_ClassesSize.end(); ++itmap)
473  {
474  os << indent << itmap->first << ": " << itmap->second << "\n";
475  }
476  }
477 
478  os << "\n" << indent << "* Training set:\n";
479  if (m_ClassesProbTraining.empty())
480  {
481  os << indent << "Not computed\n";
482  }
483  else
484  {
485  os << indent << "** Selection probability:\n";
486  for (std::map<ClassLabelType, double>::const_iterator itmap = m_ClassesProbTraining.begin();
487  itmap != m_ClassesProbTraining.end(); ++itmap)
488  {
489  os << indent << itmap->first << ": " << itmap->second << "\n";
490  }
491  os << indent << "** Number of selected samples:\n";
492  for (std::map<ClassLabelType, int>::const_iterator itmap = m_ClassesSamplesNumberTraining.begin();
493  itmap != m_ClassesSamplesNumberTraining.end(); ++itmap)
494  {
495  os << indent << itmap->first << ": " << itmap->second << "\n";
496  }
497  }
498 
499  os << "\n" << indent << "* Validation set:\n";
500  if (m_ClassesProbValidation.empty())
501  {
502  os << indent << "Not computed\n";
503  }
504  else
505  {
506  os << indent << "** Selection probability:\n";
507  for (std::map<ClassLabelType, double>::const_iterator itmap = m_ClassesProbValidation.begin();
508  itmap != m_ClassesProbValidation.end(); ++itmap)
509  {
510  os << indent << itmap->first << ": " << itmap->second << "\n";
511  }
512  os << indent << "** Number of selected samples:\n";
513  for (std::map<ClassLabelType, int>::const_iterator itmap = m_ClassesSamplesNumberValidation.begin();
514  itmap != m_ClassesSamplesNumberValidation.end(); ++itmap)
515  {
516  os << indent << itmap->first << ": " << itmap->second << "\n";
517  }
518  }
519 }
520 
521 }
522 
523 #endif
itk::Size< Monteverdi_DIMENSION > SizeType
Definition: mvdTypes.h:137
ListSampleType * GetValidationListSample()
double GetPolygonAreaInPixelsUnits(DataNodeType *polygonDataNode, ImageType *image)
const VectorDataType * GetInputVectorData() const
ListSampleType * GetTrainingListSample()
ImageType::Pointer ImagePointerType
VectorDataType::Pointer VectorDataPointerType
DataNodeType::PolygonListPointerType PolygonListPointerType
void SetInputVectorData(const VectorDataType *)
VectorDataType::DataNodeType DataNodeType
virtual const ValueType & Get() const
DataNodeType::PolygonPointerType PolygonPointerType
static ITK_CONSTEXPR_FUNC T max(const T &)
DataObjectPointerArray::size_type DataObjectPointerArraySizeType
void PrintSelf(std::ostream &os, itk::Indent indent) const override
DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) override
ImageType::RegionType TransformPhysicalRegionToIndexRegion(const RemoteSensingRegionType &region, const ImageType *image)
const ImageType * GetInput() const
virtual void SetNthOutput(DataObjectPointerArraySizeType num, DataObject *output)
ListLabelType * GetValidationListLabel()
void Fill(SizeValueType value)
void GenerateInputRequestedRegion(void) override
void SetInput(const ImageType *)
DataObject * GetOutput(const DataObjectIdentifierType &key)