otbImageFittingPolygonListFilter.hxx

Go to the documentation of this file.

/*
 * Copyright (C) 2005-2022 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 otbImageFittingPolygonListFilter_hxx
#define otbImageFittingPolygonListFilter_hxx
 
#include "otbImageFittingPolygonListFilter.h"
#include "otbPolyLineImageConstIterator.h"
#include "itkImageRegionConstIteratorWithIndex.h"
 
namespace otb
{
template <class TPath, class TImage>
ImageFittingPolygonListFilter<TPath, TImage>::ImageFittingPolygonListFilter()
{
  this->SetNumberOfRequiredInputs(2);
  this->SetNumberOfRequiredInputs(2);
  m_Radius             = 1;
  m_NumberOfIterations = 1;
}
 
template <class TPath, class TImage>
void ImageFittingPolygonListFilter<TPath, TImage>::SetInputImage(const ImageType* image)
{
  this->itk::ProcessObject::SetNthInput(1, const_cast<ImageType*>(image));
}
 
template <class TPath, class TImage>
const typename ImageFittingPolygonListFilter<TPath, TImage>::ImageType* ImageFittingPolygonListFilter<TPath, TImage>::GetInputImage(void)
{
  if (this->GetNumberOfInputs() < 1)
  {
    return nullptr;
  }
  return static_cast<const ImageType*>(this->itk::ProcessObject::GetInput(1));
}
 
// FIXME
// There is an issue here with integer and continuous indexes
// maybe we should use the itk::LineConstIterator
template <class TPath, class TImage>
void ImageFittingPolygonListFilter<TPath, TImage>::GenerateData()
{
  // I/O wiring
  ImageConstPointerType inputImagePtr = this->GetInputImage();
  const PathListType*   inputPtr      = this->GetInput();
  PathListType*         outputPtr     = this->GetOutput();
 
  typename ImageType::RegionType regionLargest = inputImagePtr->GetLargestPossibleRegion();
 
  typedef itk::ImageRegionConstIteratorWithIndex<ImageType> NeighborhoodIteratorType;
 
  typename ImageType::SizeType size;
  size[0] = 2 * m_Radius + 1;
  size[1] = 2 * m_Radius + 1;
  typename ImageType::RegionType region;
  region.SetSize(size);
  typename ImageType::IndexType start;
 
  // go through all the polygons in the list
  IteratorType it = inputPtr->Begin();
  while (it != inputPtr->End())
  {
    PathPointerType polygon = it.Get();
    if (polygon->GetVertexList()->Size() > 2)
    {
      for (unsigned int iteration = 0; iteration < m_NumberOfIterations; ++iteration)
      {
        PathPointerType             newPolygon = PathType::New();
        VertexListConstIteratorType vertexIt   = polygon->GetVertexList()->Begin();
        // We are now going to go through all the vertex, we won't start to process
        // first as we need to know the last one for that.
        VertexType firstPoint    = vertexIt.Value();
        VertexType previousPoint = vertexIt.Value();
        ++vertexIt;
        VertexType currentPoint = vertexIt.Value();
        ++vertexIt;
        while (vertexIt != polygon->GetVertexList()->End())
        {
          VertexType nextPoint = vertexIt.Value();
          {
 
            start[0] = static_cast<long int>(currentPoint[0] - m_Radius);
            start[1] = static_cast<long int>(currentPoint[1] - m_Radius);
            region.SetIndex(start);
            region.Crop(inputImagePtr->GetLargestPossibleRegion());
 
            NeighborhoodIteratorType nIt(inputImagePtr, region);
            double                   maxValue = 0.0;
            VertexType               maxPoint = currentPoint;
            nIt.GoToBegin();
            while (!nIt.IsAtEnd())
            {
              if (regionLargest.IsInside(nIt.GetIndex()))
              {
                VertexType middlePoint  = static_cast<VertexType>(nIt.GetIndex());
                double     currentValue = computeValue(inputImagePtr, middlePoint, previousPoint, nextPoint);
                if (currentValue > maxValue)
                {
                  maxValue = currentValue;
                  maxPoint = middlePoint;
                }
              }
              ++nIt;
            }
            currentPoint = maxPoint;
            newPolygon->AddVertex(maxPoint);
          }
          ++vertexIt;
          previousPoint = currentPoint;
          currentPoint  = nextPoint;
        }
        // We now need to process the last and the first point
 
        VertexType nextPoint = firstPoint;
        {
          start[0] = static_cast<long int>(currentPoint[0] - m_Radius);
          start[1] = static_cast<long int>(currentPoint[1] - m_Radius);
          region.SetIndex(start);
 
          NeighborhoodIteratorType nIt(inputImagePtr, region);
          double                   maxValue = 0.0;
          VertexType               maxPoint = currentPoint;
          nIt.GoToBegin();
 
          while (!nIt.IsAtEnd())
          {
            if (regionLargest.IsInside(nIt.GetIndex()))
            {
              VertexType middlePoint  = static_cast<VertexType>(nIt.GetIndex());
              double     currentValue = computeValue(inputImagePtr, middlePoint, previousPoint, nextPoint);
              if (currentValue > maxValue)
              {
                maxValue = currentValue;
                maxPoint = middlePoint;
              }
            }
            ++nIt;
          }
          currentPoint = maxPoint;
          newPolygon->AddVertex(maxPoint);
        }
        previousPoint = currentPoint;
        currentPoint  = firstPoint;
        vertexIt      = newPolygon->GetVertexList()->Begin();
        nextPoint     = vertexIt.Value();
 
        {
 
          start[0] = static_cast<long int>(currentPoint[0] - m_Radius);
          start[1] = static_cast<long int>(currentPoint[1] - m_Radius);
          region.SetIndex(start);
 
          NeighborhoodIteratorType nIt(inputImagePtr, region);
          double                   maxValue = 0.0;
          VertexType               maxPoint = currentPoint;
          nIt.GoToBegin();
          while (!nIt.IsAtEnd())
          {
            if (regionLargest.IsInside(nIt.GetIndex()))
            {
              VertexType middlePoint  = static_cast<VertexType>(nIt.GetIndex());
              double     currentValue = computeValue(inputImagePtr, middlePoint, previousPoint, nextPoint);
              if (currentValue > maxValue)
              {
                maxValue = currentValue;
                maxPoint = middlePoint;
              }
            }
            ++nIt;
          }
          currentPoint = maxPoint;
          newPolygon->AddVertex(maxPoint);
        }
        polygon = newPolygon; // prepare the next iteration
      }
    }
 
    outputPtr->PushBack(polygon);
 
    ++it;
  } // going through the polygon list
}
 
template <class TPath, class TImage>
double ImageFittingPolygonListFilter<TPath, TImage>::computeValue(ImageConstPointerType image, VertexType middlePoint, VertexType previousPoint,
                                                                  VertexType nextPoint) const
{
  typedef typename ImageType::IndexType IndexType;
  IndexType                             middleIndex;
  IndexType                             previousIndex;
  IndexType                             nextIndex;
  middleIndex[0]            = static_cast<long int>(middlePoint[0]);
  middleIndex[1]            = static_cast<long int>(middlePoint[1]);
  previousIndex[0]          = static_cast<long int>(previousPoint[0]);
  previousIndex[1]          = static_cast<long int>(previousPoint[1]);
  nextIndex[0]              = static_cast<long int>(nextPoint[0]);
  nextIndex[1]              = static_cast<long int>(nextPoint[1]);
  double       currentValue = 0.0;
  unsigned int count        = 0;
  // compute for first segment
  LineConstIteratorType itLineFirst(image, previousIndex, middleIndex);
  while (!itLineFirst.IsAtEnd())
  {
    currentValue += itLineFirst.Get();
    ++count;
    ++itLineFirst;
  }
 
  // compute for second segment
  LineConstIteratorType itLineSecond(image, nextIndex, middleIndex);
  while (!itLineSecond.IsAtEnd())
  {
    currentValue += itLineSecond.Get();
    ++count;
    ++itLineSecond;
  }
 
  return currentValue / count;
}
 
template <class TPath, class TImage>
void ImageFittingPolygonListFilter<TPath, TImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
{
  Superclass::PrintSelf(os, indent);
}
 
} // End namespace otb
#endif