Skip to content
Snippets Groups Projects

Resolve "Migrate code to OTB 8.x"

Merged Resolve "Migrate code to OTB 8.x"
5-migrate-code-to-otb-8-x into master
Merged Luc Hermitte requested to merge 5-migrate-code-to-otb-8-x into master
Viewing commit 3d892f9c
Prev Next
Show latest version
1 file
+ 15
− 1
Compare changes
  • Side-by-side
  • Inline
include/otbSARDEMProjectionImageFilter2.hxx 0 → 100644
+ 436
− 0
/*
* Copyright (C) 2005-2023 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.
*/
// This file is a fork of DiapOTB SARDEMProjection adapted to OTB 8 source code
// for the needs of S1Tiling
#ifndef otbSARDEMProjectionImageFilter2_hxx
#define otbSARDEMProjectionImageFilter2_hxx
#include "otbSARDEMProjectionImageFilter2.h"
#include "itkImageScanlineConstIterator.h"
#include "itkImageScanlineIterator.h"
#include "itkProgressReporter.h"
#include "itkNumericTraitsPointPixel.h"
#if OTB_VERSION_MAJOR >= 8
# include "otbGeocentricTransform.h"
#else
# if defined(__GNUC__) || defined(__clang__)
# pragma GCC diagnostic push
# pragma GCC diagnostic ignored "-Wunused-parameter"
# include "ossim/base/ossimFilename.h"
# include "ossim/base/ossimGpt.h"
# pragma GCC diagnostic pop
# else
# include "ossim/base/ossimFilename.h"
# include "ossim/base/ossimGpt.h"
# endif
#endif
#include "otbConfigurationManager.h"
#include "otbDEMHandler.h"
#include <memory>
#include <boost/optional.hpp>
#include <mutex>
#include <thread>
#include <cassert>
namespace otb
{
/**
* Constructor with default initialization
*/
template <class TImageIn, class TImageOut>
SARDEMProjectionImageFilter2<TImageIn ,TImageOut>::SARDEMProjectionImageFilter2(
std::string const& geoid_filename)
{
#if OTB_VERSION_MAJOR >= 8
auto& demHandler = DEMHandler::GetInstance();
assert(demHandler.GetGeoidFile() == geoid_filename);
#else
auto& demHandler = *DEMHandler::Instance();
if (! geoid_filename.empty())
{
const std::string isEmg96 = "egm96";
std::size_t found = geoid_filename.find(isEmg96);
if (found!=std::string::npos)
{
m_geoidEmg96 = std::make_unique<ossimGeoidEgm96>(ossimFilename(geoid_filename));
}
}
#endif
}
/**
* Set Sar Image keyWordList
*/
#if OTB_VERSION_MAJOR >= 8
template<class TImageIn, class TImageOut>
void
SARDEMProjectionImageFilter2< TImageIn ,TImageOut >
::SetSARImageMetadata(ImageMetadata sarImageMetadata)
{
m_SarImageMetadata = sarImageMetadata;
}
#else
template<class TImageIn, class TImageOut>
void
SARDEMProjectionImageFilter2< TImageIn ,TImageOut >
::SetSARImageKeyWorList(ImageKeywordlist sarImageKWL)
{
m_SarImageKwl = sarImageKWL;
}
#endif
/**
* Method GenerateOutputInformaton()
*/
template<class TImageIn, class TImageOut>
void
SARDEMProjectionImageFilter2< TImageIn, TImageOut >
::GenerateOutputInformation()
{
// Call superclass implementation
Superclass :: GenerateOutputInformation();
// Get pointers to the input and output
ImageInConstPointer inputPtr = this->GetInput();
ImageOutPointer outputPtr = this->GetOutput();
// Image KeyWordList
#if OTB_VERSION_MAJOR >= 8
auto inputKwl = inputPtr->GetImageMetadata();
#else
ImageKeywordlist inputKwl = inputPtr->GetImageKeywordlist();
#endif
// Set the number of Components for the Output VectorImage
// At Least 4 Components :
// _ C : Colunm of DEM into SAR Image
// _ L : Line of DEM into SAR Image
// _ Z : Coordonates Z
// _ Y : Coordonates Y
// (Z,Y) coordinates, defined into ossimSarSensorMode (if otb version < 8) or SarSensorMode as follows:
// Let n = |sensorPos|,
// ps2 = scalar_product(sensorPos,worldPoint)
// d = distance(sensorPos,worldPoint)
//
// Z = n - ps2/n
// Y = sqrt(d*d - Z*Z)
//
// m_withH = true => 1 component added H : high
// m_withXYZ = true => 3 Three components added Xcart Ycart Zcart : cartesien coordonates
// m_withSatPos = true => 3 Three components added : Satellite Position into cartesian
#if OTB_VERSION_MAJOR >= 8
inputKwl.Add(std::string("band.L"), std::to_string(1));
inputKwl.Add(std::string("band.C"), std::to_string(0));
inputKwl.Add(std::string("band.Z"), std::to_string(2));
inputKwl.Add(std::string("band.Y"), std::to_string(3));
#else
// Fill KeyWordList
inputKwl.AddKey("support_data.band.L", "1");
inputKwl.AddKey("support_data.band.C", "0");
inputKwl.AddKey("support_data.band.Z", "2");
inputKwl.AddKey("support_data.band.Y", "3");
#endif
m_NbComponents = 4;
if (m_withXYZ)
{
m_NbComponents += 3;
#if OTB_VERSION_MAJOR >= 8
inputKwl.Add("band.XCart", "4");
inputKwl.Add("band.YCart", "5");
inputKwl.Add("band.ZCart", "6");
#else
inputKwl.AddKey("support_data.band.XCart", "4");
inputKwl.AddKey("support_data.band.YCart", "5");
inputKwl.AddKey("support_data.band.ZCart", "6");
#endif
}
if (m_withH)
{
m_indH = m_NbComponents;
m_NbComponents += 1;
#if OTB_VERSION_MAJOR >= 8
inputKwl.Add("band.H", std::to_string(m_indH));
#else
inputKwl.AddKey("support_data.band.H", std::to_string(m_indH));
#endif
}
if (m_withSatPos)
{
m_indSatPos = m_NbComponents;
m_NbComponents += 3;
#if OTB_VERSION_MAJOR >= 8
inputKwl.Add("band.XSatPos", std::to_string(m_indSatPos));
inputKwl.Add("band.YSatPos", std::to_string(m_indSatPos+1));
inputKwl.Add("band.ZSatPos", std::to_string(m_indSatPos+2));
#else
inputKwl.AddKey("support_data.band.XSatPos", std::to_string(m_indSatPos));
inputKwl.AddKey("support_data.band.YSatPos", std::to_string(m_indSatPos+1));
inputKwl.AddKey("support_data.band.ZSatPos", std::to_string(m_indSatPos+2));
#endif
}
// Add gain and direction into kwl
#if OTB_VERSION_MAJOR >= 8
inputKwl.Add("band.DirectionToScanDEMC", std::to_string(m_directionToScanDEMC));
inputKwl.Add("band.DirectionToScanDEML", std::to_string(m_directionToScanDEML));
inputKwl.Add("band.Gain", std::to_string(m_gain));
#else
inputKwl.AddKey("support_data.band.DirectionToScanDEMC", std::to_string(m_directionToScanDEMC));
inputKwl.AddKey("support_data.band.DirectionToScanDEML", std::to_string(m_directionToScanDEML));
inputKwl.AddKey("support_data.band.Gain", std::to_string(m_gain));
#endif
outputPtr->SetNumberOfComponentsPerPixel(m_NbComponents);
#if OTB_VERSION_MAJOR >= 8
// Create and Initilaze the SarSensorModel
m_SarSensorModel = std::make_unique<SarSensorModel>(m_SarImageMetadata);
#else
// Create and Initilaze the SarSensorModelAdapter
m_SarSensorModel = SarSensorModelAdapter::New();
bool loadOk = m_SarSensorModel->LoadState(m_SarImageKwl);
if(!loadOk || !m_SarSensorModel->IsValidSensorModel())
{
itkExceptionMacro(<<"SAR image does not contain a valid SAR sensor model.");
}
#endif
// Compute the output image size, spacing and origin (same as the input)
const ImageInSizeType & inputSize = inputPtr->GetLargestPossibleRegion().GetSize();
ImageInPointType origin = inputPtr->GetOrigin();
const ImageInSpacingType & inSP = inputPtr->GetSpacing();
ImageOutRegionType outputLargestPossibleRegion = inputPtr->GetLargestPossibleRegion();
outputLargestPossibleRegion.SetSize(static_cast<ImageOutSizeType>(inputSize));
outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
outputPtr->SetOrigin(static_cast<ImageOutPointType>(origin));
outputPtr->SetSpacing(static_cast<ImageOutSpacingType>(inSP));
// Set Image KeyWordList
#if OTB_VERSION_MAJOR >= 8
outputPtr->SetImageMetadata(inputKwl);
#else
outputPtr->SetImageKeywordList(inputKwl);
#endif
}
/**
* Method GenerateInputRequestedRegion
*/
template<class TImageIn, class TImageOut>
void
SARDEMProjectionImageFilter2< TImageIn, TImageOut >
::GenerateInputRequestedRegion()
{
ImageOutRegionType outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
ImageInRegionType inputRequestedRegion = static_cast<ImageInRegionType>(outputRequestedRegion);
ImageInPointer inputPtr = const_cast< ImageInType * >( this->GetInput() );
inputPtr->SetRequestedRegion(inputRequestedRegion);
}
/**
* Method ThreadedGenerateData
*/
template<class TImageIn, class TImageOut >
void
SARDEMProjectionImageFilter2< TImageIn, TImageOut >
::ThreadedGenerateData(
ImageOutRegionType const & outputRegionForThread,
itk::ThreadIdType threadId)
{
// Compute corresponding input region (same as output here)
ImageInRegionType inputRegionForThread = static_cast<ImageInRegionType>(outputRegionForThread);
// Define/declare an iterator that will walk the input region for this
// thread.
typedef itk::ImageScanlineConstIterator< ImageInType > InputIterator;
InputIterator InIt(this->GetInput(), inputRegionForThread);
// Define/declare an iterator that will walk the output region for this
// thread. NO VECTOR IMAGE INTO ITERATOR TEMPLATE
typedef itk::ImageScanlineIterator< ImageOutType > OutputIterator;
OutputIterator OutIt(this->GetOutput(), outputRegionForThread);
// Support progress methods/callbacks
itk::ProgressReporter progress( this, threadId, outputRegionForThread.GetNumberOfPixels() );
// Transform all points
InIt.GoToBegin();
OutIt.GoToBegin();
Point3DType demGeoPoint;
Point2DType demLatLonPoint;
Point2DType col_row(0);
Point2DType y_z(0);
Point3DType xyzCart;
Point3DType satPos;
Point3DType satVel;
assert(m_NbComponents ==
4 + (m_withXYZ ? 3 : 0) + (m_withH ? 1 : 0) + (m_withSatPos ? 3 : 0));
ImageOutPixelType pixelOutput;
pixelOutput.Reserve(m_NbComponents);
#if OTB_VERSION_MAJOR >= 8
auto& demHandler = DEMHandler::GetInstance();
// demHandler.ClearElevationParameters(); // NO!!!!!
auto& demHandlerTLS = demHandler.GetHandlerForCurrentThread();
#else
auto& demHandler = *DEMHandler::Instance();
#endif
// For each Line in MNT
while ( !InIt.IsAtEnd())
{
InIt.GoToBeginOfLine();
OutIt.GoToBeginOfLine();
// For each column in MNT
while (!InIt.IsAtEndOfLine())
{
// Check if NODATA
if (InIt.Get() != m_NoData)
{
// Trnasform index to Lat/Lon Point
this->GetInput()->TransformIndexToPhysicalPoint(InIt.GetIndex(), demLatLonPoint);
demGeoPoint[0] = demLatLonPoint[0];
demGeoPoint[1] = demLatLonPoint[1];
// Get elevation from earth geoid thanks to DEMHandler or ossimgeoidEmg96
double h; // Elevation from earth geoid
#if OTB_VERSION_MAJOR >= 8
// TODO: GetGeoidHeight has very bad performances (*): GDAL
// isn't meant to be used for one coordinate at the time.
// Instead: have a multi-pass algorithm
// 1. get in a SINGLE call the geoid height for all points in
// the buffer
// 2. Compute the rest
//
// (*)
// - 96.3% of execution time spent in ThreadedGenerateData
// - 61.3% in GetGeoidHeight
// - 31.3% GDALRasterBand::RasterIO
// - 11.7% in mutex loc+unlock
// - 29.2% OCRCoordinateTransform
// - 20.2% in get_pid
//
h = GetGeoidHeight(demHandlerTLS, demLatLonPoint);
#else
if (m_geoidEmg96)
{
ossimGpt gptPt;
gptPt.lon = demLatLonPoint[0];
gptPt.lat = demLatLonPoint[1];
h = m_geoidEmg96->offsetFromEllipsoid(gptPt);
}
else
{
h = demHandler.GetHeightAboveEllipsoid(demGeoPoint[0],demGeoPoint[1]);
}
#endif
// Correct the height with earth geoid
demGeoPoint[2] = InIt.Get() + h;
// Call the method of sarSensorModel
m_SarSensorModel->WorldToLineSampleYZ(demGeoPoint, col_row, y_z);
// Fill the four channels for output
pixelOutput[0] = col_row[0];
pixelOutput[1] = col_row[1];
pixelOutput[2] = y_z[1];
pixelOutput[3] = y_z[0];
// Add channels if required
if (m_withXYZ)
{
assert(pixelOutput.Size() >= 6);
#if OTB_VERSION_MAJOR >= 8
xyzCart = Projection::WorldToEcef(demGeoPoint);
#else
SarSensorModelAdapter::WorldToCartesian(demGeoPoint, xyzCart);
#endif
pixelOutput[4] = xyzCart[0];
pixelOutput[5] = xyzCart[1];
pixelOutput[6] = xyzCart[2];
}
if (m_withH)
{
assert(pixelOutput.Size() > m_indH);
pixelOutput[m_indH] = demGeoPoint[2];
}
if (m_withSatPos)
{
assert(pixelOutput.Size() >= m_indSatPos + 3);
m_SarSensorModel->WorldToSatPositionAndVelocity(demGeoPoint, satPos, satVel);
pixelOutput[m_indSatPos] = satPos[0];
pixelOutput[m_indSatPos+1] = satPos[1];
pixelOutput[m_indSatPos+2] = satPos[2];
}
// Set the output (without iterator because no VectorImage into iterator template)
OutIt.Set(pixelOutput);
progress.CompletedPixel();
}
else
{
// Fill all channels with NoData
for (int id = 0; id < m_NbComponents; id++)
{
pixelOutput[id] = m_NoData;
}
// Set the output (without itertor because no VectorImage into iterator template)
OutIt.Set(pixelOutput);
progress.CompletedPixel();
}
// Next Col
++InIt;
++OutIt;
}
// Next Line
InIt.NextLine();
OutIt.NextLine();
}
}
} /*namespace otb*/
#endif