21 #ifndef otbVectorDataExtractROI_hxx
22 #define otbVectorDataExtractROI_hxx
27 #include "itkIdentityTransform.h"
33 #include "itkProgressReporter.h"
42 template <
class TVectorData>
50 template <
class TVectorData>
53 Superclass::PrintSelf(os, indent);
59 template <
class TVectorData>
62 this->AllocateOutputs();
63 typename VectorDataType::ConstPointer inputPtr = this->GetInput();
64 typename VectorDataType::Pointer outputPtr = this->GetOutput();
68 if (!inputPtr->GetProjectionRef().empty())
69 outputPtr->SetProjectionRef(inputPtr->GetProjectionRef());
72 this->CompareInputAndRegionProjection();
75 if (m_ProjectionNeeded)
77 otbMsgDevMacro(<<
"Reprojecting region in vector data projection");
78 this->ProjectRegionToInputVectorProjection();
82 otbMsgDevMacro(<<
"Region and vector data projection are similar");
91 typename VectorDataType::DataTreePointerType tree = outputPtr->GetDataTree();
98 newDataNode->SetNodeType(inputRoot->Get()->GetNodeType());
99 newDataNode->SetNodeId(inputRoot->Get()->GetNodeId());
101 typename InternalTreeNodeType::Pointer outputRoot = InternalTreeNodeType::New();
102 outputRoot->Set(newDataNode);
103 tree->SetRoot(outputRoot);
109 ProcessNode(inputRoot, outputRoot);
114 template <
class TVectorData>
120 for (
typename ChildrenListType::iterator it = children.begin(); it != children.end(); ++it)
122 typename InternalTreeNodeType::Pointer newContainer;
126 newDataNode->SetNodeType(dataNode->GetNodeType());
127 newDataNode->SetNodeId(dataNode->GetNodeId());
128 newDataNode->SetMetaDataDictionary(dataNode->GetMetaDataDictionary());
130 switch (dataNode->GetNodeType())
134 newContainer = InternalTreeNodeType::New();
135 newContainer->Set(newDataNode);
136 destination->AddChild(newContainer);
137 ProcessNode((*it), newContainer);
143 newContainer = InternalTreeNodeType::New();
144 newContainer->Set(newDataNode);
145 destination->AddChild(newContainer);
147 ProcessNode((*it), newContainer);
152 newContainer = InternalTreeNodeType::New();
153 newContainer->Set(newDataNode);
154 destination->AddChild(newContainer);
156 ProcessNode((*it), newContainer);
161 if (m_GeoROI.IsInside(this->PointToContinuousIndex(dataNode->GetPoint())))
163 newDataNode->SetPoint(dataNode->GetPoint());
164 newContainer = InternalTreeNodeType::New();
165 newContainer->Set(newDataNode);
166 destination->AddChild(newContainer);
173 if (this->IsLineIntersectionNotNull(dataNode->GetLine()))
175 newDataNode->SetLine(dataNode->GetLine());
176 newContainer = InternalTreeNodeType::New();
177 newContainer->Set(newDataNode);
178 destination->AddChild(newContainer);
185 if (this->IsPolygonIntersectionNotNull(dataNode->GetPolygonExteriorRing()))
187 newDataNode->SetPolygonExteriorRing(dataNode->GetPolygonExteriorRing());
188 newDataNode->SetPolygonInteriorRings(dataNode->GetPolygonInteriorRings());
189 newContainer = InternalTreeNodeType::New();
190 newContainer->Set(newDataNode);
191 destination->AddChild(newContainer);
198 newContainer = InternalTreeNodeType::New();
199 newContainer->Set(newDataNode);
200 destination->AddChild(newContainer);
202 ProcessNode((*it), newContainer);
207 newContainer = InternalTreeNodeType::New();
208 newContainer->Set(newDataNode);
209 destination->AddChild(newContainer);
211 ProcessNode((*it), newContainer);
216 newContainer = InternalTreeNodeType::New();
217 newContainer->Set(newDataNode);
218 destination->AddChild(newContainer);
220 ProcessNode((*it), newContainer);
225 newContainer = InternalTreeNodeType::New();
226 newContainer->Set(newDataNode);
227 destination->AddChild(newContainer);
229 ProcessNode((*it), newContainer);
239 template <
class TVectorData>
245 for (
unsigned int i = 0; i < polygon->GetVertexList()->Size() - 2; ++i)
250 firstVertex = polygon->GetVertexList()->GetElement(i);
251 secondVertex = polygon->GetVertexList()->GetElement(i + 1);
255 typename LineType::Pointer line = LineType::New();
256 line->AddVertex(firstVertex);
257 line->AddVertex(secondVertex);
259 if (this->IsLineIntersectionNotNull(line))
268 template <
class TVectorData>
271 RegionType lineRegion(line->GetBoundingRegion());
276 if (!lineRegion.
Crop(m_GeoROI))
283 for (
unsigned int i = 0; i < line->GetVertexList()->Size() - 1; ++i)
288 firstVertex = line->GetVertexList()->GetElement(i);
289 secondVertex = line->GetVertexList()->GetElement(i + 1);
294 firstPoint[0] = firstVertex[0];
295 firstPoint[1] = firstVertex[1];
297 secondPoint[0] = secondVertex[0];
298 secondPoint[1] = secondVertex[1];
300 if (m_GeoROI.IsInside(this->PointToContinuousIndex(firstPoint)) || m_GeoROI.IsInside(this->PointToContinuousIndex(secondPoint)))
308 if (!m_GeoROI.IsInside(this->PointToContinuousIndex(firstPoint)) && !m_GeoROI.IsInside(this->PointToContinuousIndex(secondPoint)))
311 typename LineType::Pointer tempLine = LineType::New();
312 tempLine->AddVertex(firstVertex);
313 tempLine->AddVertex(secondVertex);
316 RegionType region(tempLine->GetBoundingRegion());
317 if (region.
Crop(m_GeoROI))
332 template <
class TVectorData>
336 PointType vertexA, vertexB, vertexC, vertexD;
338 vertexA = segmentLineAB->GetVertexList()->GetElement(0);
339 vertexB = segmentLineAB->GetVertexList()->GetElement(1);
340 vertexC = segmentLineCD->GetVertexList()->GetElement(0);
341 vertexD = segmentLineCD->GetVertexList()->GetElement(1);
343 int CounterClockWiseValueWithC = CounterClockWise(vertexA, vertexB, vertexC);
344 int CounterClockWiseValueWithD = CounterClockWise(vertexA, vertexB, vertexD);
346 if (CounterClockWiseValueWithC == CounterClockWiseValueWithD)
350 int CounterClockWiseValueWithA = CounterClockWise(vertexC, vertexD, vertexA);
351 int CounterClockWiseValueWithB = CounterClockWise(vertexC, vertexD, vertexB);
353 if (CounterClockWiseValueWithA == CounterClockWiseValueWithB)
361 template <
class TVectorData>
367 SecondMinusFirstPoint = secondPoint - firstPoint;
368 ThirdMinusFirstPoint = thirdPoint - firstPoint;
370 double dX1dY2MinusdY1dX2;
371 dX1dY2MinusdY1dX2 =
static_cast<double>(SecondMinusFirstPoint[0] * ThirdMinusFirstPoint[1] - SecondMinusFirstPoint[1] * ThirdMinusFirstPoint[0]);
372 if (dX1dY2MinusdY1dX2 > 0.0)
376 if (dX1dY2MinusdY1dX2 < 0.0)
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))
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]);
396 if ((dX1dX1 + dY1dY1) < (dX2dX2 + dY2dY2))
408 template <
class TVectorData>
411 std::string regionProjection = m_ROI.GetRegionProjection();
412 std::string inputVectorProjection = this->GetInput()->GetProjectionRef();
417 if (regionProjection == inputVectorProjection)
418 m_ProjectionNeeded =
false;
420 m_ProjectionNeeded =
true;
426 template <
class TVectorData>
431 typename GenericRSTransformType::Pointer genericTransform = GenericRSTransformType::New();
435 genericTransform->SetInputProjectionRef(m_ROI.GetRegionProjection());
436 genericTransform->SetInputImageMetadata(&(m_ROI.GetImageMetadata()));
437 genericTransform->SetOutputProjectionRef(this->GetInput()->GetProjectionRef());
440 genericTransform->SetOutputOrigin(this->GetInput()->GetOrigin());
441 genericTransform->SetOutputSpacing(this->GetInput()->GetSpacing());
442 genericTransform->InstantiateTransform();
447 typename VertexListType::Pointer regionCorners = VertexListType::New();
451 point1[0] = m_ROI.GetOrigin()[0];
452 point1[1] = m_ROI.GetOrigin()[1];
455 point2[0] = m_ROI.GetOrigin()[0] + m_ROI.GetSize()[0];
456 point2[1] = m_ROI.GetOrigin()[1];
458 point3[0] = m_ROI.GetOrigin()[0] + m_ROI.GetSize()[0];
459 point3[1] = m_ROI.GetOrigin()[1] + m_ROI.GetSize()[1];
461 point4[0] = m_ROI.GetOrigin()[0];
462 point4[1] = m_ROI.GetOrigin()[1] + m_ROI.GetSize()[1];
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)));
472 m_GeoROI = this->ComputeVertexListBoundingRegion(regionCorners.GetPointer());
473 m_GeoROI.SetRegionProjection(this->GetInput()->GetProjectionRef());
479 template <
class TVectorData>
485 vertex[0] = point[0];
486 vertex[1] = point[1];
494 template <
class TVectorData>
498 double x = 0., y = 0.;
507 typename VertexListType::ConstIterator it = vertexlist->Begin();
509 if (vertexlist->Size() > 0)
511 x =
static_cast<double>(it.Value()[0]);
512 y =
static_cast<double>(it.Value()[1]);
519 while (it != vertexlist->End())
521 x =
static_cast<double>(it.Value()[0]);
522 y =
static_cast<double>(it.Value()[1]);
546 size[0] = maxId[0] - index[0];
547 size[1] = maxId[1] - index[1];