otbGCPsToRPCSensorModelImageFilter.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 otbGCPsToRPCSensorModelImageFilter_hxx
#define otbGCPsToRPCSensorModelImageFilter_hxx

#include "otbGCPsToRPCSensorModelImageFilter.h"

#include "otbRPCSolver.h"
#include "otbGenericRSTransform.h"

#include "itkMetaDataObject.h"
#include "otbMetaDataKey.h"

namespace otb
{

template <class TImage>
GCPsToRPCSensorModelImageFilter<TImage>::GCPsToRPCSensorModelImageFilter()
  : m_UseImageGCPs(false),
    m_RMSGroundError(0.),
    m_ErrorsContainer(),
    m_MeanError(0.),
    m_UseDEM(false),
    m_MeanElevation(0.),
    m_GCPsContainer(),
    m_ModelUpToDate(false)
{
  // This filter does not modify the image buffer, but only its
  // metadata.Therefore, it can be run inplace to reduce memory print.
  // CastImageFilter has InPlaceOff by default (see UnaryFunctorImgeFilter constructor)
  this->InPlaceOn();

  // Clear the GCPs container
  this->ClearGCPs();
}

template <class TImage>
GCPsToRPCSensorModelImageFilter<TImage>::~GCPsToRPCSensorModelImageFilter()
{
  // Clear the GCPs container
  this->ClearGCPs();
}


template <class TImage>
void GCPsToRPCSensorModelImageFilter<TImage>::Modified() const
{
  // Call superclass implementation
  this->Superclass::Modified();

  // Deactivate up-to-date flag
  m_ModelUpToDate = false;
}


template <class TImage>
typename GCPsToRPCSensorModelImageFilter<TImage>::GCPsContainerType& GCPsToRPCSensorModelImageFilter<TImage>::GetGCPsContainer()
{
  // return the GCPs container
  return m_GCPsContainer;
}

template <class TImage>
typename GCPsToRPCSensorModelImageFilter<TImage>::ErrorsContainerType& GCPsToRPCSensorModelImageFilter<TImage>::GetErrorsContainer()
{
  // return the GCPs container
  return m_ErrorsContainer;
}

template <class TImage>
void GCPsToRPCSensorModelImageFilter<TImage>::SetGCPsContainer(const GCPsContainerType& container)
{
  // Set the GCPs container
  m_GCPsContainer = container;
  // Call modified method
  this->Modified();
}

template <class TImage>
void GCPsToRPCSensorModelImageFilter<TImage>::AddGCP(const Point2DType& sensorPoint, const Point3DType& groundPoint)
{
  // Create a new GCP
  GCPType newGCP(sensorPoint, groundPoint);
  // Add it to the container
  m_GCPsContainer.push_back(newGCP);
  // Call the modified method */
  this->Modified();
}

template <class TImage>
void GCPsToRPCSensorModelImageFilter<TImage>::AddGCP(const Point2DType& sensorPoint, const Point2DType& groundPoint)
{
  // Create the ground point with elevation
  Point3DType groundPointWithElevation;
  // Fill it with ground point data
  groundPointWithElevation[0] = groundPoint[0];
  groundPointWithElevation[1] = groundPoint[1];

  // Check whether we are using a DEM or not
  if (m_UseDEM)
  {
    // If so, use it to get the elevation
    double height = DEMHandler::GetInstance().GetHeightAboveEllipsoid(groundPoint);
    // To avoid nan value
    if (height != height)
      height                    = 0;
    groundPointWithElevation[2] = height;
  }
  else
  {
    // Use the MeanElevation value
    groundPointWithElevation[2] = m_MeanElevation;
  }

  // Call the 3D version of the method
  this->AddGCP(sensorPoint, groundPointWithElevation);
}

template <class TImage>
void GCPsToRPCSensorModelImageFilter<TImage>::RemoveGCP(unsigned int id)
{
  if (id < m_GCPsContainer.size())
  {
    m_GCPsContainer.erase(m_GCPsContainer.begin() + id);
  }

  this->Modified();
}

template <class TImage>
void GCPsToRPCSensorModelImageFilter<TImage>::ClearGCPs()
{
  m_GCPsContainer.clear();
  this->Modified();
}

template <class TImage>
void GCPsToRPCSensorModelImageFilter<TImage>::SetUseImageGCPs(bool use)
{
  // Set the internal value
  m_UseImageGCPs = use;

  this->LoadImageGCPs();
}

template <class TImage>
void GCPsToRPCSensorModelImageFilter<TImage>::LoadImageGCPs()
{
  // First, retrieve the image pointer
  typename TImage::Pointer imagePtr = this->GetOutput();

  // Iterate on the image GCPs
  for (unsigned int i = 0; i < imagePtr->GetGCPCount(); ++i)
  {
    // Store row/col in an index
    typename TImage::IndexType index;
    index[0] = static_cast<long int>(imagePtr->GetGCPCol(i));
    index[1] = static_cast<long int>(imagePtr->GetGCPRow(i));

    // Transform to physical point
    Point2DType sensorPoint;
    imagePtr->TransformIndexToPhysicalPoint(index, sensorPoint);

    // Sensor and Ground points
    Point3DType groundPoint;

    groundPoint[0] = imagePtr->GetGCPX(i); // Lon
    groundPoint[1] = imagePtr->GetGCPY(i); // Lat
    groundPoint[2] = imagePtr->GetGCPZ(i);

    // Search image GCPs and remove them from container
    bool found = false;

    unsigned int currentGCP = 0;

    while ((currentGCP < m_GCPsContainer.size()) && !found)
    {
      if ((m_GCPsContainer[currentGCP].first == sensorPoint) && (m_GCPsContainer[currentGCP].second == groundPoint))
      {
        // If we don't use image GCPs, erase the GCP
        if (!m_UseImageGCPs)
          m_GCPsContainer.erase(m_GCPsContainer.begin() + currentGCP);

        found = true;
      }

      currentGCP++;
    }

    // If we use image GCPs, add it to the container
    if (m_UseImageGCPs && !found)
      this->AddGCP(sensorPoint, groundPoint);
  }

  this->Modified();
}

template <class TImage>
void GCPsToRPCSensorModelImageFilter<TImage>::ComputeErrors()
{
  using RSTransformType = GenericRSTransform<double, 2, 2>;

  RSTransformType::Pointer rsTransform = RSTransformType::New();
  rsTransform->SetInputImageMetadata(&m_ImageMetadata);

  rsTransform->InstantiateTransform();

  double sum  = 0.;
  m_MeanError = 0.;

  // Clear Error container
  m_ErrorsContainer.clear();

  for (unsigned int i = 0; i < m_GCPsContainer.size(); ++i)
  {
    // GCP value
    const auto & sensorPoint = m_GCPsContainer[i].first;
    const auto & groundPoint = m_GCPsContainer[i].second;

    // Compute Transform
    Point3DType groundPointTemp, sensorPointTemp;
    sensorPointTemp[0] = sensorPoint[0];
    sensorPointTemp[1] = sensorPoint[1];

    auto outPoint = rsTransform->TransformPoint(sensorPoint);

    Point2DType groundPoint2D;

    groundPoint2D[0] = groundPoint[0];
    groundPoint2D[1] = groundPoint[1];

    double error = outPoint.EuclideanDistanceTo(outPoint);

    // Add error to the container
    m_ErrorsContainer.push_back(error);

    // Compute mean error
    sum += error;
  }

  m_MeanError = sum / static_cast<double>(m_ErrorsContainer.size());
}

template <class TImage>
void GCPsToRPCSensorModelImageFilter<TImage>::GenerateOutputInformation()
{
  // First, call the superclass implementation
  Superclass::GenerateOutputInformation();

  // First, retrieve the image pointer
  auto imagePtr = this->GetOutput();

  if (!m_ModelUpToDate)
  {
    double rmsError;

    Projection::RPCParam rpcParams;
    otb::RPCSolver::Solve(m_GCPsContainer, rmsError, rpcParams);

    // Retrieve the residual ground error
    m_RMSGroundError = rmsError;

    m_ImageMetadata = imagePtr->GetImageMetadata();
    m_ImageMetadata.Add(MDGeom::RPC, rpcParams);

    // Compute errors
    this->ComputeErrors();

    m_ModelUpToDate = true;
  }

  imagePtr->SetImageMetadata(m_ImageMetadata);
}

template <class TImage>
void GCPsToRPCSensorModelImageFilter<TImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
{
  Superclass::PrintSelf(os, indent);
  os << indent << "UseImageGCPs: " << (m_UseImageGCPs ? "yes" : "no") << std::endl;
  os << indent << "UseDEM: " << (m_UseDEM ? "yes" : "no") << std::endl;
  if (!m_UseDEM)
  {
    os << indent << "MeanElevation: " << m_MeanElevation << std::endl;
  }
  os << indent << "RMS ground error: " << m_RMSGroundError << std::endl;
}

} // end of namespace otb

#endif