otbVectorDataExtractROI.hxx

Go to the documentation of this file.

/*
 * Copyright (C) 2005-2024 Centre National d'Etudes Spatiales (CNES)
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
 
#ifndef otbVectorDataExtractROI_hxx
#define otbVectorDataExtractROI_hxx
 
#include "otbVectorDataExtractROI.h"
 
#include "otbGenericMapProjection.h"
#include "itkIdentityTransform.h"
#include "otbGenericRSTransform.h"
 
#include "otbObjectList.h"
#include "otbMacro.h"
 
#include "itkProgressReporter.h"
#include "otbStopwatch.h"
 
namespace otb
{
 
template <class TVectorData>
VectorDataExtractROI<TVectorData>::VectorDataExtractROI() : m_ProjectionNeeded(false), m_ROI(), m_Kept(0)
{
}
 
template <class TVectorData>
void VectorDataExtractROI<TVectorData>::PrintSelf(std::ostream& os, itk::Indent indent) const
{
  Superclass::PrintSelf(os, indent);
}
 
template <class TVectorData>
void VectorDataExtractROI<TVectorData>::GenerateData(void)
{
  this->AllocateOutputs();
  typename VectorDataType::ConstPointer inputPtr  = this->GetInput();
  typename VectorDataType::Pointer      outputPtr = this->GetOutput();
 
  // Find out the projection needed
  if (!inputPtr->GetProjectionRef().empty())
    outputPtr->SetProjectionRef(inputPtr->GetProjectionRef());
 
  this->CompareInputAndRegionProjection();
 
  if (m_ProjectionNeeded)
  {
    otbMsgDevMacro(<< "Reprojecting region in vector data projection");
    this->ProjectRegionToInputVectorProjection();
  }
  else
  {
    otbMsgDevMacro(<< "Region and vector data projection are similar");
    m_GeoROI = m_ROI;
  }
 
  otbMsgDevMacro(<< "ROI: " << this->m_ROI);
  otbMsgDevMacro(<< "GeoROI: " << this->m_GeoROI);
 
  // Create the output tree root
  // DataNodePointerType outputRoot = DataNodeType::New();
  // outputRoot->SetNodeType(inputPtr->GetRoot()->GetNodeType());
  // outputRoot->SetNodeId(inputPtr->GetRoot()->GetNodeId());
  // outputPtr->SetRoot(outputRoot);
 
  m_Kept = 0;
 
  // Start recursive processing
  otb::Stopwatch chrono = otb::Stopwatch::StartNew();
  ProcessNode(inputPtr,inputPtr->GetRoot(),outputPtr,outputPtr->GetRoot());
  chrono.Stop();
  otbMsgDevMacro(<< "VectorDataExtractROI: " << m_Kept << " features processed in " << chrono.GetElapsedMilliseconds() << " ms.");
} /*End GenerateData()*/
 
template <class TVectorData>
void VectorDataExtractROI<TVectorData>::ProcessNode(VectorDataConstPointerType inputVdata, DataNodePointerType source, VectorDataPointerType outputVdata, DataNodePointerType destination)
{
  // Get the children list from the input node
  ChildrenListType children = inputVdata->GetChildrenList(source);
  // For each child
  for (typename ChildrenListType::iterator it = children.begin(); it != children.end(); ++it)
  {
    DataNodePointerType dataNode    = (*it);
    DataNodePointerType newDataNode = DataNodeType::New();
    newDataNode->SetNodeType(dataNode->GetNodeType());
    newDataNode->SetNodeId(dataNode->GetNodeId());
    newDataNode->SetMetaDataDictionary(dataNode->GetMetaDataDictionary());
 
    switch (dataNode->GetNodeType())
    {
    case ROOT:
    {
      //outputVdata->Add(newDataNode,destination);
      //ProcessNode(inputVdata,(*it), outputVdata, newDataNode);
      //++m_Kept;
      break;
    }
    case DOCUMENT:
    {
      outputVdata->Add(newDataNode,destination);
      ++m_Kept;
      ProcessNode(inputVdata,(*it), outputVdata, newDataNode);
      break;
    }
    case FOLDER:
    {
      outputVdata->Add(newDataNode,destination);
      ++m_Kept;
      ProcessNode(inputVdata,(*it), outputVdata, newDataNode);
      break;
    }
    case FEATURE_POINT:
    {
      if (m_GeoROI.IsInside(this->PointToContinuousIndex(dataNode->GetPoint())))
      {
        newDataNode->SetPoint(dataNode->GetPoint());
        outputVdata->Add(newDataNode,destination);
        ++m_Kept;
      }
      break;
    }
    case FEATURE_LINE:
    {
      if (this->IsLineIntersectionNotNull(dataNode->GetLine()))
      {
        newDataNode->SetLine(dataNode->GetLine());
        outputVdata->Add(newDataNode,destination);
        ++m_Kept;
      }
      break;
    }
    case FEATURE_POLYGON:
    {
      if (this->IsPolygonIntersectionNotNull(dataNode->GetPolygonExteriorRing()))
      {
        newDataNode->SetPolygonExteriorRing(dataNode->GetPolygonExteriorRing());
        newDataNode->SetPolygonInteriorRings(dataNode->GetPolygonInteriorRings());
        outputVdata->Add(newDataNode,destination);
        ++m_Kept;
      }
      break;
    }
    case FEATURE_MULTIPOINT:
    {
      outputVdata->Add(newDataNode,destination);
      ++m_Kept;
      ProcessNode(inputVdata,(*it), outputVdata, newDataNode);
      break;
    }
    case FEATURE_MULTILINE:
    {
      outputVdata->Add(newDataNode,destination);
      ++m_Kept;
      ProcessNode(inputVdata,(*it), outputVdata, newDataNode);
      break;
    }
    case FEATURE_MULTIPOLYGON:
    {
      outputVdata->Add(newDataNode,destination);
      ++m_Kept;
      ProcessNode(inputVdata,(*it), outputVdata, newDataNode);
      break;
    }
    case FEATURE_COLLECTION:
    {
      outputVdata->Add(newDataNode,destination);
      ++m_Kept;
      ProcessNode(inputVdata,(*it), outputVdata, newDataNode);
      break;
    }
    }
  }
}
 
template <class TVectorData>
bool VectorDataExtractROI<TVectorData>::IsPolygonIntersectionNotNull(PolygonPointerType polygon)
{
  // Get the VertexList
  // -2 cause we don't want the last point
  // which is the same as the first one (closed Ring)
  for (unsigned int i = 0; i < polygon->GetVertexList()->Size() - 2; ++i)
  {
    // Get the components of the polygon 2 by 2
    VertexType firstVertex;
    VertexType secondVertex;
    firstVertex  = polygon->GetVertexList()->GetElement(i);
    secondVertex = polygon->GetVertexList()->GetElement(i + 1);
 
    // Build a line with each two vertex
    typename LineType::Pointer line = LineType::New();
    line->AddVertex(firstVertex);
    line->AddVertex(secondVertex);
 
    if (this->IsLineIntersectionNotNull(line))
      return true;
  }
  return false;
}
 
template <class TVectorData>
bool VectorDataExtractROI<TVectorData>::IsLineIntersectionNotNull(LinePointerType line)
{
  RegionType lineRegion(line->GetBoundingRegion());
 
  // if the line bounding box have a null
  // intersection with the geoROI
  // the line and the region do not intersect
  if (!lineRegion.Crop(m_GeoROI))
  {
    return false;
  }
  else
  {
    // Get the VertexList
    for (unsigned int i = 0; i < line->GetVertexList()->Size() - 1; ++i)
    {
      // Get the components of the line 2 by 2
      VertexType firstVertex;
      VertexType secondVertex;
      firstVertex  = line->GetVertexList()->GetElement(i);
      secondVertex = line->GetVertexList()->GetElement(i + 1);
 
      //  -------------------
      // Case 1 : Check if one of the two points are in the region
      PointType firstPoint, secondPoint;
      firstPoint[0] = firstVertex[0];
      firstPoint[1] = firstVertex[1];
 
      secondPoint[0] = secondVertex[0];
      secondPoint[1] = secondVertex[1];
 
      if (m_GeoROI.IsInside(this->PointToContinuousIndex(firstPoint)) || m_GeoROI.IsInside(this->PointToContinuousIndex(secondPoint)))
      {
        return true;
      }
      else
      {
        //  -------------------
        // Case 2 : If any of the point is in the region
        if (!m_GeoROI.IsInside(this->PointToContinuousIndex(firstPoint)) && !m_GeoROI.IsInside(this->PointToContinuousIndex(secondPoint)))
        {
          // Build a line with each two vertex
          typename LineType::Pointer tempLine = LineType::New();
          tempLine->AddVertex(firstVertex);
          tempLine->AddVertex(secondVertex);
 
          // Check if the intersection is not null
          RegionType region(tempLine->GetBoundingRegion());
          if (region.Crop(m_GeoROI))
            return true;
 
          //  -------------------
          // TODO : check if the segment cut
          //        one of the region edges
        }
      }
    }
  }
 
  return false;
}
 
 
template <class TVectorData>
bool VectorDataExtractROI<TVectorData>::IsSegmentIntersectSegment(LinePointerType segmentLineAB, LinePointerType segmentLineCD)
{
 
  PointType vertexA, vertexB, vertexC, vertexD;
 
  vertexA = segmentLineAB->GetVertexList()->GetElement(0);
  vertexB = segmentLineAB->GetVertexList()->GetElement(1);
  vertexC = segmentLineCD->GetVertexList()->GetElement(0);
  vertexD = segmentLineCD->GetVertexList()->GetElement(1);
 
  int CounterClockWiseValueWithC = CounterClockWise(vertexA, vertexB, vertexC);
  int CounterClockWiseValueWithD = CounterClockWise(vertexA, vertexB, vertexD);
 
  if (CounterClockWiseValueWithC == CounterClockWiseValueWithD)
  {
    return false;
  }
  int CounterClockWiseValueWithA = CounterClockWise(vertexC, vertexD, vertexA);
  int CounterClockWiseValueWithB = CounterClockWise(vertexC, vertexD, vertexB);
 
  if (CounterClockWiseValueWithA == CounterClockWiseValueWithB)
  {
    return false;
  }
 
  return true;
}
 
template <class TVectorData>
int VectorDataExtractROI<TVectorData>::CounterClockWise(PointType firstPoint, PointType secondPoint, PointType thirdPoint)
{
  PointType SecondMinusFirstPoint;
  PointType ThirdMinusFirstPoint;
 
  SecondMinusFirstPoint = secondPoint - firstPoint;
  ThirdMinusFirstPoint  = thirdPoint - firstPoint;
 
  double dX1dY2MinusdY1dX2;
  dX1dY2MinusdY1dX2 = static_cast<double>(SecondMinusFirstPoint[0] * ThirdMinusFirstPoint[1] - SecondMinusFirstPoint[1] * ThirdMinusFirstPoint[0]);
  if (dX1dY2MinusdY1dX2 > 0.0)
  {
    return 1;
  }
  if (dX1dY2MinusdY1dX2 < 0.0)
  {
    return -1;
  }
 
  double dX1dX2;
  double dY1dY2;
  dX1dX2 = static_cast<double>(SecondMinusFirstPoint[0] * ThirdMinusFirstPoint[0]);
  dY1dY2 = static_cast<double>(SecondMinusFirstPoint[1] * ThirdMinusFirstPoint[1]);
  if ((dX1dX2 < 0.0) || (dY1dY2 < 0.0))
  {
    return -1;
  }
 
  double dX1dX1, dX2dX2, dY1dY1, dY2dY2;
  dX1dX1 = static_cast<double>(SecondMinusFirstPoint[0] * SecondMinusFirstPoint[0]);
  dX2dX2 = static_cast<double>(ThirdMinusFirstPoint[0] * ThirdMinusFirstPoint[0]);
  dY1dY1 = static_cast<double>(SecondMinusFirstPoint[1] * SecondMinusFirstPoint[1]);
  dY2dY2 = static_cast<double>(ThirdMinusFirstPoint[1] * ThirdMinusFirstPoint[1]);
 
  if ((dX1dX1 + dY1dY1) < (dX2dX2 + dY2dY2))
  {
    return 1;
  }
 
  return 0;
}
 
 
template <class TVectorData>
void VectorDataExtractROI<TVectorData>::CompareInputAndRegionProjection()
{
  std::string regionProjection      = m_ROI.GetRegionProjection();
  std::string inputVectorProjection = this->GetInput()->GetProjectionRef();
 
  // FIXME: the string comparison is not sufficient to say that two
  // projections are different
  if (regionProjection == inputVectorProjection)
    m_ProjectionNeeded = false;
  else
    m_ProjectionNeeded = true;
}
 
template <class TVectorData>
void VectorDataExtractROI<TVectorData>::ProjectRegionToInputVectorProjection()
{
  /* Use the RS Generic projection */
  typedef otb::GenericRSTransform<>        GenericRSTransformType;
  typename GenericRSTransformType::Pointer genericTransform = GenericRSTransformType::New();
 
  genericTransform->SetInputProjectionRef(m_ROI.GetRegionProjection());
  genericTransform->SetInputImageMetadata(&(m_ROI.GetImageMetadata()));
  genericTransform->SetOutputProjectionRef(this->GetInput()->GetProjectionRef());
  //TODO: const itk::MetaDataDictionary& inputDict = this->GetInput()->GetMetaDataDictionary();
  //TODO: genericTransform->SetOutputImageMetadata(this->GetInput()->GetImageMetadata());
  genericTransform->SetOutputOrigin(this->GetInput()->GetOrigin());
  genericTransform->SetOutputSpacing(this->GetInput()->GetSpacing());
  genericTransform->InstantiateTransform();
 
  otbMsgDevMacro(<< genericTransform);
 
  typename VertexListType::Pointer regionCorners = VertexListType::New();
  ProjPointType                    point1, point2, point3, point4;
 
  point1[0] = m_ROI.GetOrigin()[0];
  point1[1] = m_ROI.GetOrigin()[1];
 
  point2[0] = m_ROI.GetOrigin()[0] + m_ROI.GetSize()[0];
  point2[1] = m_ROI.GetOrigin()[1];
 
  point3[0] = m_ROI.GetOrigin()[0] + m_ROI.GetSize()[0];
  point3[1] = m_ROI.GetOrigin()[1] + m_ROI.GetSize()[1];
 
  point4[0] = m_ROI.GetOrigin()[0];
  point4[1] = m_ROI.GetOrigin()[1] + m_ROI.GetSize()[1];
 
  regionCorners->InsertElement(regionCorners->Size(), this->PointToContinuousIndex(genericTransform->TransformPoint(point1)));
  regionCorners->InsertElement(regionCorners->Size(), this->PointToContinuousIndex(genericTransform->TransformPoint(point2)));
  regionCorners->InsertElement(regionCorners->Size(), this->PointToContinuousIndex(genericTransform->TransformPoint(point3)));
  regionCorners->InsertElement(regionCorners->Size(), this->PointToContinuousIndex(genericTransform->TransformPoint(point4)));
 
  m_GeoROI = this->ComputeVertexListBoundingRegion(regionCorners.GetPointer());
  m_GeoROI.SetRegionProjection(this->GetInput()->GetProjectionRef());
}
 
template <class TVectorData>
typename VectorDataExtractROI<TVectorData>::VertexType VectorDataExtractROI<TVectorData>::PointToContinuousIndex(ProjPointType point)
{
 
  VertexType vertex;
 
  vertex[0] = point[0];
  vertex[1] = point[1];
 
  return vertex;
}
 
template <class TVectorData>
typename VectorDataExtractROI<TVectorData>::RegionType
VectorDataExtractROI<TVectorData>::ComputeVertexListBoundingRegion(typename VertexListType::ConstPointer vertexlist)
{
  double    x = 0., y = 0.;
  IndexType index;
  IndexType maxId;
  SizeType  size;
 
  index.Fill(0.);
  maxId.Fill(0.);
  size.Fill(0.);
 
  typename VertexListType::ConstIterator it = vertexlist->Begin();
 
  if (vertexlist->Size() > 0)
  {
    x        = static_cast<double>(it.Value()[0]);
    y        = static_cast<double>(it.Value()[1]);
    index[0] = x;
    index[1] = y;
    maxId[0] = x;
    maxId[1] = y;
 
    ++it;
    while (it != vertexlist->End())
    {
      x = static_cast<double>(it.Value()[0]);
      y = static_cast<double>(it.Value()[1]);
 
      // Index search
      if (x < index[0])
      {
        index[0] = x;
      }
      if (y > index[1])
      {
        index[1] = y;
      }
      // Max Id search for size computation
      if (x > maxId[0])
      {
        maxId[0] = x;
      }
      if (y < maxId[1])
      {
        maxId[1] = y;
      }
 
      ++it;
    }
 
    size[0] = maxId[0] - index[0];
    size[1] = maxId[1] - index[1];
  }
 
  RegionType region;
  region.SetSize(size);
  region.SetOrigin(index);
 
  return region;
}
 
} // end namespace otb
 
#endif