OTB  7.2.0
Orfeo Toolbox
otbVectorDataExtractROI.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2020 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 
73 
76  {
77  otbMsgDevMacro(<< "Reprojecting region in vector data projection");
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->SetInputKeywordList(m_ROI.GetKeywordList());
437  genericTransform->SetOutputProjectionRef(this->GetInput()->GetProjectionRef());
438  const itk::MetaDataDictionary& inputDict = this->GetInput()->GetMetaDataDictionary();
439  genericTransform->SetOutputDictionary(inputDict);
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
bool IsSegmentIntersectSegment(LinePointerType segmentLineAB, LinePointerType segmentLineCD)
RegionType::IndexType IndexType
virtual OutputVectorDataType * GetOutput(void)
bool IsInside(const IndexType &index) const
DataNodeType::PolygonPointerType PolygonPointerType
InternalTreeNodeType::ChildrenListType ChildrenListType
static Stopwatch StartNew()
DataNodeType::PointType PointType
VectorDataType::DataTreeType::TreeNodeType InternalTreeNodeType
PolygonType::VertexType VertexType
itk::Point< typename VertexType::CoordRepType, IndexType::IndexDimension > ProjPointType
Stopwatch timer.
Definition: otbStopwatch.h:41
DurationType GetElapsedMilliseconds() const
DataNodeType::Pointer DataNodePointerType
virtual void ProjectRegionToInputVectorProjection()
virtual VertexType PointToContinuousIndex(ProjPointType point)
const IndexType & GetOrigin() const
virtual RegionType ComputeVertexListBoundingRegion(typename VertexListType::ConstPointer vertexlist)
void SetOrigin(const IndexType &index)
void PrintSelf(std::ostream &os, itk::Indent indent) const override
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
virtual bool IsPolygonIntersectionNotNull(PolygonPointerType polygon)
This is the class to handle generic remote sensing transform.
int CounterClockWise(PointType firstPoint, PointType secondPoint, PointType thirdPoint)
void SetRegionProjection(const std::string &projection)
virtual void ProcessNode(InternalTreeNodeType *source, InternalTreeNodeType *destination)
DataNodeType::LinePointerType LinePointerType
virtual bool IsLineIntersectionNotNull(LinePointerType line)
void SetSize(const SizeType &size)
const ImageKeywordlist & GetKeywordList() const
void PrintSelf(std::ostream &os, itk::Indent indent) const override
const SizeType & GetSize() const
#define otbMsgDevMacro(x)
Definition: otbMacro.h:64