otbImageToGenericRSOutputParameters.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 otbImageToGenericRSOutputParameters_hxx
#define otbImageToGenericRSOutputParameters_hxx
 
 
#include "otbImageToGenericRSOutputParameters.h"
#include "itkMacro.h"
#include "itkContinuousIndex.h"
 
namespace otb
{
 
template <class TImage>
ImageToGenericRSOutputParameters<TImage>::ImageToGenericRSOutputParameters()
{
  m_Transform                = GenericRSTransformType::New();
  m_ForceSpacing             = false;
  m_ForceSize                = false;
  m_EstimateIsotropicSpacing = false;
}
 
template <class TImage>
void ImageToGenericRSOutputParameters<TImage>::Compute()
{
  // Do some checks : exceptions if Null or empty projectionRef and empty keywordlist
  if (m_Input.IsNull())
    itkExceptionMacro(<< "The input is null , please set a non null input image");
 
  if (!m_Input->GetImageMetadata().HasSensorGeometry() && !m_Input->GetImageMetadata().HasProjectedGeometry())
    itkExceptionMacro(<< "No information in the metadata, please set an image with non empty metadata");
 
  // First Call to UpdateTransform : Initialize with the input image
  // information
  this->UpdateTransform();
 
  // Estimate the Output image Extent
  this->EstimateOutputImageExtent();
 
  // Estimate the Output Spacing
  if (!m_ForceSpacing)
    this->EstimateOutputSpacing();
 
  // Finally Estimate the Output Size
  this->EstimateOutputSize();
 
  // Estimate the Output Origin
  this->EstimateOutputOrigin();
}
 
 
template <class TImage>
void ImageToGenericRSOutputParameters<TImage>::UpdateTransform()
{
  m_Transform->SetOutputImageMetadata(&(this->GetInput()->GetImageMetadata()));
  m_Transform->SetOutputProjectionRef(this->GetInput()->GetProjectionRef());
  m_Transform->InstantiateTransform();
}
 
 
template <class TImage>
void ImageToGenericRSOutputParameters<TImage>::EstimateOutputImageExtent()
{
  // Get the inverse transform again : used later
  GenericRSTransformPointerType invTransform = GenericRSTransformType::New();
  m_Transform->GetInverse(invTransform);
 
  // Compute the 4 corners in the cartographic coordinate system
  std::vector<itk::ContinuousIndex<double, 2>> vindex;
  std::vector<PointType> voutput;
 
  itk::ContinuousIndex<double, 2> index1(m_Input->GetLargestPossibleRegion().GetIndex());
  index1[0] += -0.5;
  index1[1] += -0.5;
  itk::ContinuousIndex<double, 2> index2(index1);
  itk::ContinuousIndex<double, 2> index3(index1);
  itk::ContinuousIndex<double, 2> index4(index1);
 
  // Image size
  SizeType size = m_Input->GetLargestPossibleRegion().GetSize();
 
  // project the 4 corners
  index2[0] += size[0];
  index3[0] += size[0];
  index3[1] += size[1];
  index4[1] += size[1];
 
  vindex.push_back(index1);
  vindex.push_back(index2);
  vindex.push_back(index3);
  vindex.push_back(index4);
 
  for (unsigned int i = 0; i < vindex.size(); ++i)
  {
    PointType physicalPoint;
    m_Input->TransformContinuousIndexToPhysicalPoint(vindex[i], physicalPoint);
    voutput.push_back(invTransform->TransformPoint(physicalPoint));
  }
 
  // Compute the boundaries
  double minX = voutput[0][0];
  double maxX = voutput[0][0];
  double minY = voutput[0][1];
  double maxY = voutput[0][1];
 
  for (unsigned int i = 0; i < voutput.size(); ++i)
  {
    // Origins
    if (minX > voutput[i][0])
      minX = voutput[i][0];
    if (minY > voutput[i][1])
      minY = voutput[i][1];
 
    // Sizes
    if (maxX < voutput[i][0])
      maxX = voutput[i][0];
 
    if (maxY < voutput[i][1])
      maxY = voutput[i][1];
  }
 
  // Edit the output image extent type
  m_OutputExtent.maxX = maxX;
  m_OutputExtent.minX = minX;
  m_OutputExtent.maxY = maxY;
  m_OutputExtent.minY = minY;
}
 
 
template <class TImage>
void ImageToGenericRSOutputParameters<TImage>::EstimateOutputOrigin()
{
  // Set the output origin in carto
  // projection
  PointType origin;
  origin[0] = m_OutputExtent.minX + 0.5 * this->GetOutputSpacing()[0];
  origin[1] = m_OutputExtent.maxY + 0.5 * this->GetOutputSpacing()[1];
  this->SetOutputOrigin(origin);
}
 
template <class TImage>
void ImageToGenericRSOutputParameters<TImage>::EstimateOutputSpacing()
{
  // Compute the output size
  double sizeCartoX = std::abs(m_OutputExtent.maxX - m_OutputExtent.minX);
  double sizeCartoY = std::abs(m_OutputExtent.minY - m_OutputExtent.maxY);
 
  PointType o, oX, oY;
  o[0] = m_OutputExtent.minX;
  o[1] = m_OutputExtent.maxY;
 
  oX = o;
  oY = o;
 
  oX[0] += sizeCartoX;
  oY[1] -= sizeCartoY;
 
  // Transform back into the input image
  PointType io  = m_Transform->TransformPoint(o);
  PointType ioX = m_Transform->TransformPoint(oX);
  PointType ioY = m_Transform->TransformPoint(oY);
 
  // Transform to indices
  IndexType ioIndex, ioXIndex, ioYIndex;
  m_Input->TransformPhysicalPointToIndex(io, ioIndex);
  m_Input->TransformPhysicalPointToIndex(ioX, ioXIndex);
  m_Input->TransformPhysicalPointToIndex(ioY, ioYIndex);
 
  // Evaluate Ox and Oy length in number of pixels
  double OxLength, OyLength;
 
  OxLength = std::sqrt(std::pow((double)ioIndex[0] - (double)ioXIndex[0], 2) + std::pow((double)ioIndex[1] - (double)ioXIndex[1], 2));
 
  OyLength = std::sqrt(std::pow((double)ioIndex[0] - (double)ioYIndex[0], 2) + std::pow((double)ioIndex[1] - (double)ioYIndex[1], 2));
 
  // Evaluate spacing
  SpacingType outputSpacing;
 
 
  if (m_EstimateIsotropicSpacing)
  {
    double isotropicSpacing = std::min(sizeCartoX / OxLength, sizeCartoY / OyLength);
    outputSpacing[0]        = isotropicSpacing;
    outputSpacing[1]        = -isotropicSpacing;
  }
  else
  {
    outputSpacing[0] = sizeCartoX / OxLength;
    outputSpacing[1] = -sizeCartoY / OyLength;
  }
 
  this->SetOutputSpacing(outputSpacing);
}
 
template <class TImage>
void ImageToGenericRSOutputParameters<TImage>::EstimateOutputSize()
{
  // Compute the output size
  double sizeCartoX = std::abs(m_OutputExtent.maxX - m_OutputExtent.minX);
  double sizeCartoY = std::abs(m_OutputExtent.minY - m_OutputExtent.maxY);
 
  // Evaluate output size
  SizeType outputSize;
  outputSize[0] = static_cast<unsigned int>(std::floor(std::abs(sizeCartoX / this->GetOutputSpacing()[0])));
  outputSize[1] = static_cast<unsigned int>(std::floor(std::abs(sizeCartoY / this->GetOutputSpacing()[1])));
 
  // if ForceSizeTo used don't update the output size with the value
  // computed : the value is computed to update the spacing knowing
  // the forced size and the computed one.
  if (!m_ForceSize)
    this->SetOutputSize(outputSize);
  else
  {
    // Compute the spacing knowing the size
    SpacingType outputSpacing;
    outputSpacing[0] = this->GetOutputSpacing()[0] * outputSize[0] / this->GetOutputSize()[0];
    outputSpacing[1] = this->GetOutputSpacing()[1] * outputSize[1] / this->GetOutputSize()[1];
    this->SetOutputSpacing(outputSpacing);
  }
  this->UpdateTransform();
}
}
 
#endif