Resolve "Migrate code to OTB 8.x"

Closes #5 (closed)

include/otbSARDEMProjectionImageFilter2.hxx

+/*
+ * 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