OTB  9.0.0
Orfeo Toolbox
otbVectorDataExtractROI.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 otbVectorDataExtractROI_hxx
22 #define otbVectorDataExtractROI_hxx
23 
25 
27 #include "itkIdentityTransform.h"
28 #include "otbGenericRSTransform.h"
29 
30 #include "otbObjectList.h"
31 #include "otbMacro.h"
32 
33 #include "itkProgressReporter.h"
34 #include "otbStopwatch.h"
35 
36 namespace otb
37 {
38 
42 template <class TVectorData>
43 VectorDataExtractROI<TVectorData>::VectorDataExtractROI() : m_ProjectionNeeded(false), m_ROI(), m_Kept(0)
44 {
45 }
46 
50 template <class TVectorData>
51 void VectorDataExtractROI<TVectorData>::PrintSelf(std::ostream& os, itk::Indent indent) const
52 {
53  Superclass::PrintSelf(os, indent);
54 }
55 
59 template <class TVectorData>
61 {
62  this->AllocateOutputs();
63  typename VectorDataType::ConstPointer inputPtr = this->GetInput();
64  typename VectorDataType::Pointer outputPtr = this->GetOutput();
66 
67  // Find out the projection needed
68  if (!inputPtr->GetProjectionRef().empty())
69  outputPtr->SetProjectionRef(inputPtr->GetProjectionRef());
70 
72  this->CompareInputAndRegionProjection();
73 
75  if (m_ProjectionNeeded)
76  {
77  otbMsgDevMacro(<< "Reprojecting region in vector data projection");
78  this->ProjectRegionToInputVectorProjection();
79  }
80  else
81  {
82  otbMsgDevMacro(<< "Region and vector data projection are similar");
83  m_GeoROI = m_ROI;
84  }
86 
87  otbMsgDevMacro(<< "ROI: " << this->m_ROI);
88  otbMsgDevMacro(<< "GeoROI: " << this->m_GeoROI);
89 
90  // Retrieve the output tree
91  typename VectorDataType::DataTreePointerType tree = outputPtr->GetDataTree();
92 
93  // Get the input tree root
94  InternalTreeNodeType* inputRoot = const_cast<InternalTreeNodeType*>(inputPtr->GetDataTree()->GetRoot());
95 
96  // Create the output tree root
97  DataNodePointerType newDataNode = DataNodeType::New();
98  newDataNode->SetNodeType(inputRoot->Get()->GetNodeType());
99  newDataNode->SetNodeId(inputRoot->Get()->GetNodeId());
100 
101  typename InternalTreeNodeType::Pointer outputRoot = InternalTreeNodeType::New();
102  outputRoot->Set(newDataNode);
103  tree->SetRoot(outputRoot);
104 
105  m_Kept = 0;
106 
107  // Start recursive processing
109  ProcessNode(inputRoot, outputRoot);
110  chrono.Stop();
111  otbMsgDevMacro(<< "VectorDataExtractROI: " << m_Kept << " features processed in " << chrono.GetElapsedMilliseconds() << " ms.");
112 } /*End GenerateData()*/
113 
114 template <class TVectorData>
116 {
117  // Get the children list from the input node
118  ChildrenListType children = source->GetChildrenList();
119  // For each child
120  for (typename ChildrenListType::iterator it = children.begin(); it != children.end(); ++it)
121  {
122  typename InternalTreeNodeType::Pointer newContainer;
123 
124  DataNodePointerType dataNode = (*it)->Get();
125  DataNodePointerType newDataNode = DataNodeType::New();
126  newDataNode->SetNodeType(dataNode->GetNodeType());
127  newDataNode->SetNodeId(dataNode->GetNodeId());
128  newDataNode->SetMetaDataDictionary(dataNode->GetMetaDataDictionary());
129 
130  switch (dataNode->GetNodeType())
131  {
132  case ROOT:
133  {
134  newContainer = InternalTreeNodeType::New();
135  newContainer->Set(newDataNode);
136  destination->AddChild(newContainer);
137  ProcessNode((*it), newContainer);
138  ++m_Kept;
139  break;
140  }
141  case DOCUMENT:
142  {
143  newContainer = InternalTreeNodeType::New();
144  newContainer->Set(newDataNode);
145  destination->AddChild(newContainer);
146  ++m_Kept;
147  ProcessNode((*it), newContainer);
148  break;
149  }
150  case FOLDER:
151  {
152  newContainer = InternalTreeNodeType::New();
153  newContainer->Set(newDataNode);
154  destination->AddChild(newContainer);
155  ++m_Kept;
156  ProcessNode((*it), newContainer);
157  break;
158  }
159  case FEATURE_POINT:
160  {
161  if (m_GeoROI.IsInside(this->PointToContinuousIndex(dataNode->GetPoint())))
162  {
163  newDataNode->SetPoint(dataNode->GetPoint());
164  newContainer = InternalTreeNodeType::New();
165  newContainer->Set(newDataNode);
166  destination->AddChild(newContainer);
167  ++m_Kept;
168  }
169  break;
170  }
171  case FEATURE_LINE:
172  {
173  if (this->IsLineIntersectionNotNull(dataNode->GetLine()))
174  {
175  newDataNode->SetLine(dataNode->GetLine());
176  newContainer = InternalTreeNodeType::New();
177  newContainer->Set(newDataNode);
178  destination->AddChild(newContainer);
179  ++m_Kept;
180  }
181  break;
182  }
183  case FEATURE_POLYGON:
184  {
185  if (this->IsPolygonIntersectionNotNull(dataNode->GetPolygonExteriorRing()))
186  {
187  newDataNode->SetPolygonExteriorRing(dataNode->GetPolygonExteriorRing());
188  newDataNode->SetPolygonInteriorRings(dataNode->GetPolygonInteriorRings());
189  newContainer = InternalTreeNodeType::New();
190  newContainer->Set(newDataNode);
191  destination->AddChild(newContainer);
192  ++m_Kept;
193  }
194  break;
195  }
196  case FEATURE_MULTIPOINT:
197  {
198  newContainer = InternalTreeNodeType::New();
199  newContainer->Set(newDataNode);
200  destination->AddChild(newContainer);
201  ++m_Kept;
202  ProcessNode((*it), newContainer);
203  break;
204  }
205  case FEATURE_MULTILINE:
206  {
207  newContainer = InternalTreeNodeType::New();
208  newContainer->Set(newDataNode);
209  destination->AddChild(newContainer);
210  ++m_Kept;
211  ProcessNode((*it), newContainer);
212  break;
213  }
215  {
216  newContainer = InternalTreeNodeType::New();
217  newContainer->Set(newDataNode);
218  destination->AddChild(newContainer);
219  ++m_Kept;
220  ProcessNode((*it), newContainer);
221  break;
222  }
223  case FEATURE_COLLECTION:
224  {
225  newContainer = InternalTreeNodeType::New();
226  newContainer->Set(newDataNode);
227  destination->AddChild(newContainer);
228  ++m_Kept;
229  ProcessNode((*it), newContainer);
230  break;
231  }
232  }
233  }
234 }
235 
239 template <class TVectorData>
241 {
242  // Get the VertexList
243  // -2 cause we don't want the last point
244  // which is the same as the first one (closed Ring)
245  for (unsigned int i = 0; i < polygon->GetVertexList()->Size() - 2; ++i)
246  {
247  // Get the components of the polygon 2 by 2
248  VertexType firstVertex;
249  VertexType secondVertex;
250  firstVertex = polygon->GetVertexList()->GetElement(i);
251  secondVertex = polygon->GetVertexList()->GetElement(i + 1);
253 
254  // Build a line with each two vertex
255  typename LineType::Pointer line = LineType::New();
256  line->AddVertex(firstVertex);
257  line->AddVertex(secondVertex);
258 
259  if (this->IsLineIntersectionNotNull(line))
260  return true;
261  }
262  return false;
263 }
264 
268 template <class TVectorData>
270 {
271  RegionType lineRegion(line->GetBoundingRegion());
272 
273  // if the line bounding box have a null
274  // intersection with the geoROI
275  // the line and the region do not intersect
276  if (!lineRegion.Crop(m_GeoROI))
277  {
278  return false;
279  }
280  else
281  {
282  // Get the VertexList
283  for (unsigned int i = 0; i < line->GetVertexList()->Size() - 1; ++i)
284  {
285  // Get the components of the line 2 by 2
286  VertexType firstVertex;
287  VertexType secondVertex;
288  firstVertex = line->GetVertexList()->GetElement(i);
289  secondVertex = line->GetVertexList()->GetElement(i + 1);
290 
291  // -------------------
292  // Case 1 : Check if one of the two points are in the region
293  PointType firstPoint, secondPoint;
294  firstPoint[0] = firstVertex[0];
295  firstPoint[1] = firstVertex[1];
296 
297  secondPoint[0] = secondVertex[0];
298  secondPoint[1] = secondVertex[1];
299 
300  if (m_GeoROI.IsInside(this->PointToContinuousIndex(firstPoint)) || m_GeoROI.IsInside(this->PointToContinuousIndex(secondPoint)))
301  {
302  return true;
303  }
304  else
305  {
306  // -------------------
307  // Case 2 : If any of the point is in the region
308  if (!m_GeoROI.IsInside(this->PointToContinuousIndex(firstPoint)) && !m_GeoROI.IsInside(this->PointToContinuousIndex(secondPoint)))
309  {
310  // Build a line with each two vertex
311  typename LineType::Pointer tempLine = LineType::New();
312  tempLine->AddVertex(firstVertex);
313  tempLine->AddVertex(secondVertex);
314 
315  // Check if the intersection is not null
316  RegionType region(tempLine->GetBoundingRegion());
317  if (region.Crop(m_GeoROI))
318  return true;
319 
320  // -------------------
321  // TODO : check if the segment cut
322  // one of the region edges
323  }
324  }
325  }
326  }
327 
328  return false;
329 }
330 
331 
332 template <class TVectorData>
334 {
335 
336  PointType vertexA, vertexB, vertexC, vertexD;
337 
338  vertexA = segmentLineAB->GetVertexList()->GetElement(0);
339  vertexB = segmentLineAB->GetVertexList()->GetElement(1);
340  vertexC = segmentLineCD->GetVertexList()->GetElement(0);
341  vertexD = segmentLineCD->GetVertexList()->GetElement(1);
342 
343  int CounterClockWiseValueWithC = CounterClockWise(vertexA, vertexB, vertexC);
344  int CounterClockWiseValueWithD = CounterClockWise(vertexA, vertexB, vertexD);
345 
346  if (CounterClockWiseValueWithC == CounterClockWiseValueWithD)
347  {
348  return false;
349  }
350  int CounterClockWiseValueWithA = CounterClockWise(vertexC, vertexD, vertexA);
351  int CounterClockWiseValueWithB = CounterClockWise(vertexC, vertexD, vertexB);
352 
353  if (CounterClockWiseValueWithA == CounterClockWiseValueWithB)
354  {
355  return false;
356  }
357 
358  return true;
359 }
360 
361 template <class TVectorData>
363 {
364  PointType SecondMinusFirstPoint;
365  PointType ThirdMinusFirstPoint;
366 
367  SecondMinusFirstPoint = secondPoint - firstPoint;
368  ThirdMinusFirstPoint = thirdPoint - firstPoint;
369 
370  double dX1dY2MinusdY1dX2;
371  dX1dY2MinusdY1dX2 = static_cast<double>(SecondMinusFirstPoint[0] * ThirdMinusFirstPoint[1] - SecondMinusFirstPoint[1] * ThirdMinusFirstPoint[0]);
372  if (dX1dY2MinusdY1dX2 > 0.0)
373  {
374  return 1;
375  }
376  if (dX1dY2MinusdY1dX2 < 0.0)
377  {
378  return -1;
379  }
380 
381  double dX1dX2;
382  double dY1dY2;
383  dX1dX2 = static_cast<double>(SecondMinusFirstPoint[0] * ThirdMinusFirstPoint[0]);
384  dY1dY2 = static_cast<double>(SecondMinusFirstPoint[1] * ThirdMinusFirstPoint[1]);
385  if ((dX1dX2 < 0.0) || (dY1dY2 < 0.0))
386  {
387  return -1;
388  }
389 
390  double dX1dX1, dX2dX2, dY1dY1, dY2dY2;
391  dX1dX1 = static_cast<double>(SecondMinusFirstPoint[0] * SecondMinusFirstPoint[0]);
392  dX2dX2 = static_cast<double>(ThirdMinusFirstPoint[0] * ThirdMinusFirstPoint[0]);
393  dY1dY1 = static_cast<double>(SecondMinusFirstPoint[1] * SecondMinusFirstPoint[1]);
394  dY2dY2 = static_cast<double>(ThirdMinusFirstPoint[1] * ThirdMinusFirstPoint[1]);
395 
396  if ((dX1dX1 + dY1dY1) < (dX2dX2 + dY2dY2))
397  {
398  return 1;
399  }
400 
401  return 0;
402 }
403 
404 
408 template <class TVectorData>
410 {
411  std::string regionProjection = m_ROI.GetRegionProjection();
412  std::string inputVectorProjection = this->GetInput()->GetProjectionRef();
414 
415  // FIXME: the string comparison is not sufficient to say that two
416  // projections are different
417  if (regionProjection == inputVectorProjection)
418  m_ProjectionNeeded = false;
419  else
420  m_ProjectionNeeded = true;
421 }
422 
426 template <class TVectorData>
428 {
429  /* Use the RS Generic projection */
430  typedef otb::GenericRSTransform<> GenericRSTransformType;
431  typename GenericRSTransformType::Pointer genericTransform = GenericRSTransformType::New();
433 
435  genericTransform->SetInputProjectionRef(m_ROI.GetRegionProjection());
436  genericTransform->SetInputImageMetadata(&(m_ROI.GetImageMetadata()));
437  genericTransform->SetOutputProjectionRef(this->GetInput()->GetProjectionRef());
438  //TODO: const itk::MetaDataDictionary& inputDict = this->GetInput()->GetMetaDataDictionary();
439  //TODO: genericTransform->SetOutputImageMetadata(this->GetInput()->GetImageMetadata());
440  genericTransform->SetOutputOrigin(this->GetInput()->GetOrigin());
441  genericTransform->SetOutputSpacing(this->GetInput()->GetSpacing());
442  genericTransform->InstantiateTransform();
444 
445  otbMsgDevMacro(<< genericTransform);
446 
447  typename VertexListType::Pointer regionCorners = VertexListType::New();
448  ProjPointType point1, point2, point3, point4;
449 
451  point1[0] = m_ROI.GetOrigin()[0];
452  point1[1] = m_ROI.GetOrigin()[1];
454 
455  point2[0] = m_ROI.GetOrigin()[0] + m_ROI.GetSize()[0];
456  point2[1] = m_ROI.GetOrigin()[1];
457 
458  point3[0] = m_ROI.GetOrigin()[0] + m_ROI.GetSize()[0];
459  point3[1] = m_ROI.GetOrigin()[1] + m_ROI.GetSize()[1];
460 
461  point4[0] = m_ROI.GetOrigin()[0];
462  point4[1] = m_ROI.GetOrigin()[1] + m_ROI.GetSize()[1];
463 
465  regionCorners->InsertElement(regionCorners->Size(), this->PointToContinuousIndex(genericTransform->TransformPoint(point1)));
466  regionCorners->InsertElement(regionCorners->Size(), this->PointToContinuousIndex(genericTransform->TransformPoint(point2)));
467  regionCorners->InsertElement(regionCorners->Size(), this->PointToContinuousIndex(genericTransform->TransformPoint(point3)));
468  regionCorners->InsertElement(regionCorners->Size(), this->PointToContinuousIndex(genericTransform->TransformPoint(point4)));
470 
472  m_GeoROI = this->ComputeVertexListBoundingRegion(regionCorners.GetPointer());
473  m_GeoROI.SetRegionProjection(this->GetInput()->GetProjectionRef());
474 }
475 
479 template <class TVectorData>
481 {
482 
483  VertexType vertex;
484 
485  vertex[0] = point[0];
486  vertex[1] = point[1];
487 
488  return vertex;
489 }
490 
494 template <class TVectorData>
496 VectorDataExtractROI<TVectorData>::ComputeVertexListBoundingRegion(typename VertexListType::ConstPointer vertexlist)
497 {
498  double x = 0., y = 0.;
499  IndexType index;
500  IndexType maxId;
501  SizeType size;
502 
503  index.Fill(0.);
504  maxId.Fill(0.);
505  size.Fill(0.);
506 
507  typename VertexListType::ConstIterator it = vertexlist->Begin();
508 
509  if (vertexlist->Size() > 0)
510  {
511  x = static_cast<double>(it.Value()[0]);
512  y = static_cast<double>(it.Value()[1]);
513  index[0] = x;
514  index[1] = y;
515  maxId[0] = x;
516  maxId[1] = y;
517 
518  ++it;
519  while (it != vertexlist->End())
520  {
521  x = static_cast<double>(it.Value()[0]);
522  y = static_cast<double>(it.Value()[1]);
523 
524  // Index search
525  if (x < index[0])
526  {
527  index[0] = x;
528  }
529  if (y > index[1])
530  {
531  index[1] = y;
532  }
533  // Max Id search for size computation
534  if (x > maxId[0])
535  {
536  maxId[0] = x;
537  }
538  if (y < maxId[1])
539  {
540  maxId[1] = y;
541  }
542 
543  ++it;
544  }
545 
546  size[0] = maxId[0] - index[0];
547  size[1] = maxId[1] - index[1];
548  }
549 
550  RegionType region;
551  region.SetSize(size);
552  region.SetOrigin(index);
553 
554  return region;
555 }
556 
557 } // end namespace otb
558 
559 #endif
otb::FEATURE_POINT
@ FEATURE_POINT
Definition: otbDataNode.h:43
otb::RemoteSensingRegion::SetOrigin
void SetOrigin(const IndexType &index)
Definition: otbRemoteSensingRegion.h:132
otb::FEATURE_MULTIPOLYGON
@ FEATURE_MULTIPOLYGON
Definition: otbDataNode.h:48
otb::DOCUMENT
@ DOCUMENT
Definition: otbDataNode.h:41
otb::VectorDataExtractROI::ProcessNode
virtual void ProcessNode(InternalTreeNodeType *source, InternalTreeNodeType *destination)
Definition: otbVectorDataExtractROI.hxx:115
otb::VectorDataExtractROI::GenerateData
void GenerateData(void) override
Definition: otbVectorDataExtractROI.hxx:60
otb::VectorDataExtractROI::PrintSelf
void PrintSelf(std::ostream &os, itk::Indent indent) const override
Definition: otbVectorDataExtractROI.hxx:51
otbGenericMapProjection.h
otb::RemoteSensingRegion::SetSize
void SetSize(const SizeType &size)
Definition: otbRemoteSensingRegion.h:164
otb
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
Definition: otbJoinContainer.h:32
otb::FEATURE_LINE
@ FEATURE_LINE
Definition: otbDataNode.h:44
otbMacro.h
otb::VectorDataExtractROI::CompareInputAndRegionProjection
virtual void CompareInputAndRegionProjection()
Definition: otbVectorDataExtractROI.hxx:409
otb::VectorDataExtractROI::IsSegmentIntersectSegment
bool IsSegmentIntersectSegment(LinePointerType segmentLineAB, LinePointerType segmentLineCD)
Definition: otbVectorDataExtractROI.hxx:333
otb::FOLDER
@ FOLDER
Definition: otbDataNode.h:42
otb::VectorDataExtractROI::InternalTreeNodeType
VectorDataType::DataTreeType::TreeNodeType InternalTreeNodeType
Definition: otbVectorDataExtractROI.h:93
otb::VectorDataExtractROI::ComputeVertexListBoundingRegion
virtual RegionType ComputeVertexListBoundingRegion(typename VertexListType::ConstPointer vertexlist)
Definition: otbVectorDataExtractROI.hxx:496
otb::VectorDataExtractROI::ProjPointType
itk::Point< typename VertexType::CoordRepType, IndexType::IndexDimension > ProjPointType
Definition: otbVectorDataExtractROI.h:90
otb::Stopwatch
Stopwatch timer.
Definition: otbStopwatch.h:41
otb::VectorDataExtractROI::SizeType
RegionType::SizeType SizeType
Definition: otbVectorDataExtractROI.h:88
otb::VectorDataExtractROI::IsLineIntersectionNotNull
virtual bool IsLineIntersectionNotNull(LinePointerType line)
Definition: otbVectorDataExtractROI.hxx:269
otb::VectorDataExtractROI::IsPolygonIntersectionNotNull
virtual bool IsPolygonIntersectionNotNull(PolygonPointerType polygon)
Definition: otbVectorDataExtractROI.hxx:240
otb::VectorDataExtractROI::DataNodePointerType
DataNodeType::Pointer DataNodePointerType
Definition: otbVectorDataExtractROI.h:72
otb::GenericRSTransform
This is the class to handle generic remote sensing transform.
Definition: otbGenericRSTransform.h:57
otb::RemoteSensingRegion::Crop
bool Crop(const Self &region)
Definition: otbRemoteSensingRegion.h:270
otb::VectorDataExtractROI::ChildrenListType
InternalTreeNodeType::ChildrenListType ChildrenListType
Definition: otbVectorDataExtractROI.h:94
otb::Stopwatch::StartNew
static Stopwatch StartNew()
otbObjectList.h
otb::RemoteSensingRegion< typename VertexType::CoordRepType >
otb::FEATURE_POLYGON
@ FEATURE_POLYGON
Definition: otbDataNode.h:45
otb::VectorDataExtractROI::IndexType
RegionType::IndexType IndexType
Definition: otbVectorDataExtractROI.h:87
otb::VectorDataExtractROI::CounterClockWise
int CounterClockWise(PointType firstPoint, PointType secondPoint, PointType thirdPoint)
Definition: otbVectorDataExtractROI.hxx:362
otb::FEATURE_MULTIPOINT
@ FEATURE_MULTIPOINT
Definition: otbDataNode.h:46
otb::VectorDataExtractROI::PointToContinuousIndex
virtual VertexType PointToContinuousIndex(ProjPointType point)
Definition: otbVectorDataExtractROI.hxx:480
otbGenericRSTransform.h
otb::VectorDataExtractROI::VectorDataExtractROI
VectorDataExtractROI()
Definition: otbVectorDataExtractROI.hxx:43
otb::Stopwatch::Stop
void Stop()
otbMsgDevMacro
#define otbMsgDevMacro(x)
Definition: otbMacro.h:64
otb::VectorDataExtractROI::ProjectRegionToInputVectorProjection
virtual void ProjectRegionToInputVectorProjection()
Definition: otbVectorDataExtractROI.hxx:427
otb::Stopwatch::GetElapsedMilliseconds
DurationType GetElapsedMilliseconds() const
otb::FEATURE_COLLECTION
@ FEATURE_COLLECTION
Definition: otbDataNode.h:49
otb::FEATURE_MULTILINE
@ FEATURE_MULTILINE
Definition: otbDataNode.h:47
otb::ROOT
@ ROOT
Definition: otbDataNode.h:40
otb::VectorDataExtractROI::PolygonPointerType
DataNodeType::PolygonPointerType PolygonPointerType
Definition: otbVectorDataExtractROI.h:74
otb::VectorDataExtractROI::VertexType
PolygonType::VertexType VertexType
Definition: otbVectorDataExtractROI.h:81
otbVectorDataExtractROI.h
otbStopwatch.h
otb::VectorDataExtractROI::LinePointerType
DataNodeType::LinePointerType LinePointerType
Definition: otbVectorDataExtractROI.h:77
otb::VectorDataExtractROI::PointType
DataNodeType::PointType PointType
Definition: otbVectorDataExtractROI.h:78