Merge branch '2-adapt-diapotb-sarcartesianmeanestimation-to-keep-shadows' into 'master'

Resolve "Adapt DiapOTB SARCartesianMeanEstimation to keep shadows"

Closes #2

See merge request !1

app/CMakeLists.txt

 #
-# Copyright (C) 2005-2020 Centre National d'Etudes Spatiales (CNES)
+# Copyright (C) 2005-2022 Centre National d'Etudes Spatiales (CNES)
 #
 # This file is part of S1Tiling remote module for Orfeo Toolbox
 #
@@ -32,3 +32,9 @@ otb_create_application(
   NAME           SARComputeLocalIncidenceAngle
   SOURCES        otbSARComputeLocalIncidenceAngle.cxx
   LINK_LIBRARIES ${otb-module} ${${otb-module}_LIBRARIES})
+
+otb_create_application(
+  NAME           SARCartesianMeanEstimation2
+  SOURCES        otbSARCartesianMeanEstimation2.cxx
+  LINK_LIBRARIES ${${otb-module}_LIBRARIES})
+

app/otbSARCartesianMeanEstimation2.cxx

+/*
+ * 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.
+ */
+
+#include "otbWrapperApplication.h"
+#include "otbWrapperApplicationFactory.h"
+
+#include "otbSARDEMPolygonsAnalysisImageFilter2.h"
+#include "otbSARCartesianMeanFunctor2.h"
+
+#include "otbWrapperOutputImageParameter.h"
+
+#include <string>
+#include <sstream>
+
+namespace otb
+{
+namespace Wrapper
+{
+
+class SARCartesianMeanEstimation2 : public Application
+{
+public:
+  using Self                     = SARCartesianMeanEstimation2;
+  using Pointer                  = itk::SmartPointer<Self>;
+
+  itkNewMacro(Self);
+  itkTypeMacro(SARCartesianMeanEstimation2, otb::Wrapper::Application);
+
+  // Pixels
+  using ImageInPixelType         = typename FloatVectorImageType::PixelType;
+  using ImageOutPixelType        = typename FloatImageType::PixelType;
+
+  // Function
+  using PolygonsFunctorType      = otb::Function::SARPolygonsFunctor<ImageInPixelType, ImageOutPixelType>;
+  using CartesianMeanFunctorType = otb::Function::SARCartesianMeanFunctor<ImageInPixelType, ImageInPixelType>;
+
+  // Filters
+  using FilterType               = otb::SARDEMPolygonsAnalysisImageFilter < FloatVectorImageType, FloatVectorImageType, FloatImageType, ComplexFloatImageType, CartesianMeanFunctorType >;
+
+private:
+  void DoInit() override
+  {
+    SetName("SARCartesianMeanEstimation2");
+    SetDescription("SAR Cartesian mean estimation thanks to the associated DEM.");
+
+    SetDocLongDescription("This application estimates a simulated cartesian mean image thanks to a DEM file.");
+
+    //Optional descriptors
+    SetDocLimitations("Only Sentinel 1 (IW and StripMap mode) and Cosmo products are supported for now.");
+    SetDocAuthors("OTB-Team");
+    SetDocSeeAlso(" ");
+    AddDocTag(Tags::SAR);
+    AddDocTag("DiapOTB");
+
+    //Parameter declarations
+    AddParameter(ParameterType_InputImage,  "indemproj",   "Input vector of DEM projected into SAR geometry");
+    SetParameterDescription("indemproj", "Input vector image for cartesian mean estimation.");
+
+    AddParameter(ParameterType_InputImage,  "indem",   "Input DEM");
+    SetParameterDescription("indem", "DEM to extract DEM geometry.");
+
+    AddParameter(ParameterType_InputImage,  "insar",   "Input SAR image");
+    SetParameterDescription("insar", "SAR Image to extract SAR geometry.");
+
+    AddParameter(ParameterType_Int,  "indirectiondemc", "Range direction for DEM scan");
+    SetParameterDescription("indirectiondemc", "Range direction for DEM scan.");
+    SetDefaultParameterInt("indirectiondemc", 1);
+    MandatoryOff("indirectiondemc");
+
+    AddParameter(ParameterType_Int,  "indirectiondeml", "Azimut direction for DEM scan");
+    SetParameterDescription("indirectiondeml", "Azimut direction for DEM scan.");
+    SetDefaultParameterInt("indirectiondeml", 1);
+    MandatoryOff("indirectiondeml");
+
+    AddParameter(ParameterType_Int, "mlran", "Averaging on distance (output geometry)");
+    SetParameterDescription("mlran","Averaging on distance (output geometry)");
+    SetDefaultParameterInt("mlran", 3);
+    SetMinimumParameterIntValue("mlran", 1);
+    MandatoryOff("mlran");
+
+    AddParameter(ParameterType_Int, "mlazi", "Averaging on azimut (output geometry)");
+    SetParameterDescription("mlazi","Averaging on azimut (output geometry)");
+    SetDefaultParameterInt("mlazi", 3);
+    SetMinimumParameterIntValue("mlazi", 1);
+    MandatoryOff("mlazi");
+
+    AddParameter(ParameterType_OutputImage, "out", "Output cartesian (mean) Image for DEM Projection");
+    SetParameterDescription("out","Output cartesian (mean) Image for DEM Projection."
+        " The first 3 bands contain X, Y and Z cartesian coordinates of the point in original SAR geometry."
+        " The last band is a mask band: 0: XYZ coordinates cannot be computed; 1: coordinates are nominally computed; 2: coordinates are estimated, but point is estimated in shadow"
+        );
+
+
+    AddRAMParameter();
+
+    SetDocExampleParameterValue("indemproj","CLZY_S21E055.tiff");
+    SetDocExampleParameterValue("indem","S21E055.hgt");
+    SetDocExampleParameterValue("insar","s1a-s4-slc-vv-20160818t014650-20160818t014715-012648-013db1-002_SLC.tiff");
+    SetDocExampleParameterValue("out","s1a-s4-simu-cartMean.tiff");
+  }
+
+  void DoUpdateParameters() override
+  {
+    // Nothing to do here : all parameters are independent
+  }
+
+  void DoExecute() override
+  {
+    // Get numeric parameters : ML factors and directions for DEM scan
+    int mlRan = GetParameterInt("mlran");
+    int mlAzi = GetParameterInt("mlazi");
+    int DEMScanDirectionC = GetParameterInt("indirectiondemc");
+    int DEMScanDirectionL = GetParameterInt("indirectiondeml");
+
+    otbAppLogINFO(<<"ML Range : "<<mlRan);
+    otbAppLogINFO(<<"ML Azimut : "<<mlAzi);
+    otbAppLogINFO(<<"Direction (DEM scan) in range : "<<DEMScanDirectionC);
+    otbAppLogINFO(<<"Direction (DEM scan) in azimut : "<<DEMScanDirectionL);
+
+    // Start the first pipelines (Read SAR and DEM image metedata)
+    ComplexFloatImageType::Pointer SARPtr = GetParameterComplexFloatImage("insar");
+    FloatImageType::Pointer DEMPtr = GetParameterFloatImage("indem");
+
+    // CartesianMeanEstimation Filter (use SARDEMPolygonsAnalysisImageFilter with SARCartesianMeanEstimationFunctor
+    // to estimate cartesian coordonates mean image)
+    FilterType::Pointer filterCartesianMeanEstimation = FilterType::New();
+    m_Ref.push_back(filterCartesianMeanEstimation.GetPointer());
+    filterCartesianMeanEstimation->SetSARImageKeyWorList(SARPtr->GetImageKeywordlist());
+    filterCartesianMeanEstimation->SetSARImagePtr(SARPtr);
+    filterCartesianMeanEstimation->SetDEMImagePtr(DEMPtr);
+    filterCartesianMeanEstimation->SetDEMInformation(0, DEMScanDirectionC, DEMScanDirectionL);
+    filterCartesianMeanEstimation->initializeMarginAndRSTransform();
+
+    // Function Type : SARCartesianMeanEstimationFunctor
+    int nbThreads = filterCartesianMeanEstimation->GetNumberOfThreads();
+    CartesianMeanFunctorType::Pointer functor =  CartesianMeanFunctorType::New(nbThreads, mlRan, mlAzi,
+        SARPtr->GetOrigin()[0]-0.5);
+    filterCartesianMeanEstimation->SetFunctorPtr(functor);
+
+
+    // Define the main pipeline (controlled with extended FileName)
+    std::string origin_FileName = GetParameterString("out");
+
+    // Check if FileName is extended (with the ? caracter)
+    // If not extended then override the FileName
+    if (origin_FileName.find("?") == std::string::npos && !origin_FileName.empty())
+    {
+      std::string extendedFileName = origin_FileName;
+
+      // Get the ram value (in MB)
+      int ram = GetParameterInt("ram");
+
+      // Define with the ram value, the number of lines for the streaming
+      int nbColSAR = SARPtr->GetLargestPossibleRegion().GetSize()[0];
+      int nbLinesSAR = SARPtr->GetLargestPossibleRegion().GetSize()[1];
+      // To determine the number of lines :
+      // nbColSAR * nbLinesStreaming * sizeof(OutputPixel) = RamValue/2 (/2 to be sure)
+      // OutputPixel are float => sizeof(OutputPixel) = 4 (Bytes)
+      long int ram_In_KBytes = (ram/2) * 1024;
+      long int nbLinesStreaming = (ram_In_KBytes / (nbColSAR)) * (1024/4);
+
+      // Check the value of nbLinesStreaming
+      int nbLinesStreamingMax = 10000;
+      if (nbLinesStreamingMax > nbLinesSAR)
+      {
+        nbLinesStreamingMax = nbLinesSAR;
+      }
+      if (nbLinesStreaming <= 0 || nbLinesStreaming >  nbLinesStreamingMax)
+      {
+        nbLinesStreaming =  nbLinesStreamingMax;
+      }
+
+      // Construct the extendedPart
+      std::ostringstream os;
+      os << "?&streaming:type=stripped&streaming:sizevalue=" << nbLinesStreaming;
+
+      // Add the extendedPart
+      std::string extendedPart = os.str();
+      extendedFileName.append(extendedPart);
+
+      // Set the new FileName with extended options
+      SetParameterString("out", extendedFileName);
+    }
+
+    // Execute the main Pipeline
+    filterCartesianMeanEstimation->SetInput(GetParameterImage("indemproj"));
+    SetParameterOutputImage("out", filterCartesianMeanEstimation->GetOutput());
+  }
+  // Vector for filters
+  std::vector<itk::ProcessObject::Pointer> m_Ref;
+};
+
+}
+}
+
+OTB_APPLICATION_EXPORT(otb::Wrapper::SARCartesianMeanEstimation2)

include/equal_range_interval.h

+/*
+ * 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 equal_range_interval_h
+#define equal_range_interval_h
+
+#include <algorithm>  // std::lower_bound, std::upper_bound
+#include <cassert>
+#include <functional> // std::less
+#include <utility>    // std::pair
+
+namespace otb
+{
+/**
+ * Adaptation of `std::equal_range` for finding interval bounds.
+ * \tparam RandomAccessIterator Iterator type that should model `std::random_access_iterator`.
+ * \tparam T                    Type of the interval bounds
+ * \tparam Comp                 Comparison function, defaults to `std::less<>`
+ * \param[in] first  Start of the range where bounds are searched
+ * \param[in] last   Past-end of the range where bounds are searched
+ * \param[in] lo     Low bound of the interval
+ * \param[in] hi     High bound of the interval
+ * \param[in] cmp    Comparison predicate
+ *
+ * \return a pair of iterators `[it_low, it_high[` such as:
+ * - `∀ it ∈ [first, it_low[  =>  cmp(*it, lo)`
+ * - `∀ it ∈ [it_low, last[   => !cmp(*it, lo)`
+ * - `∀ it ∈ [first, it_high[ => !cmp(hi, *it)`
+ * - `∀ it ∈ [it_high, last[  =>  cmp(hi, *it)`
+ *
+ * Somehow, `it_low` is computed with `std::lower_bound`, and `it_hi` is
+ * computed with `std::upper_bound`.
+ *
+ * \pre `cmp(lo, hi)`
+ * \pre `[first, last[` shall be partitioned according to `cmp` predicate
+ * around `lo` and around `high`. If the range is sorted according to `cmp`,
+ * then this also implies it's partitionned.
+ * \throw None (it depends on `Cmp`...)
+ */
+template <typename RandomAccessIterator, typename T,
+         typename Comp = std::less<typename RandomAccessIterator::value_type>>
+inline
+std::pair<RandomAccessIterator, RandomAccessIterator>
+equal_range_interval(
+    RandomAccessIterator first, RandomAccessIterator last,
+    T lo, T hi,
+    Comp cmp = Comp{})
+{
+  // assert(std::is_sorted(first, last, cmp));
+  auto len = std::distance(first, last);
+  while (len > 0) {
+    auto half = len / 2;
+    auto middle = first + half;
+    if (cmp(*middle, lo)) {
+      first = middle + 1;
+      len -= half + 1;
+    } else if (cmp(hi, *middle)) {
+      len = half;
+    } else {
+      auto left  = std::lower_bound(first, middle, lo, cmp);
+      auto right = std::upper_bound(middle+1, first+len, hi, cmp);
+      return {left, right};
+    }
+  }
+  return {first, first};
+}
+} // namespace otb
+
+#endif  // equal_range_interval_h

include/otbPositionHelpers.h

@@ -23,6 +23,7 @@
 
 #include "otbRepack.h"
 #include "otbSarSensorModelAdapter.h"
+#include <itkContinuousIndex.h>
 
 namespace otb
 {
@@ -36,6 +37,19 @@ auto TransformIndexToPhysicalPoint(ImageType & im, typename ImageType::IndexType
   return p;
 }
 
+template <typename IndexRep, typename ImageType>
+inline
+auto TransformPhysicalPointToContinuousIndex(ImageType & im, typename ImageType::PointType const& p)
+{
+  // There is no ImageType::ContinuousIndex in otb images
+  // => work around it
+  // static_assert(std::is_same<typename ImageType::IndexValueType, double>::value, "");
+  static_assert(ImageType::ImageDimension == 2, "");
+  itk::ContinuousIndex<IndexRep, ImageType::ImageDimension> index;
+  im.TransformPhysicalPointToContinuousIndex(p, index);
+  return index;
+}
+
 template <typename Point3DType>
 inline
 auto WorldToCartesian(Point3DType const& lat_lon_hgt)

include/otbSARCartesianMeanFunctor2.h

+/*
+ * 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 otbSARCartesianMeanFunctor2_h
+#define otbSARCartesianMeanFunctor2_h
+
+#include "otbSARPolygonsFunctor2.h"
+
+#include "itkObject.h"
+#include "itkObjectFactory.h"
+#include <cassert>
+
+namespace otb
+{
+namespace Function
+{
+/** \class SARCartesianMeanFunctor
+ * \brief Catesian mean estimation.
+ *
+ * This functor is an implementation of SARPolygonsFunctor. It estimates the Cartesian coordonates mean image.
+ *
+ * The characteristics are the following :
+ * _ Output Geometry : ML (defined by ML factors)
+ * _ Seven required components : L, C, Y, Z, XCart, YCart and ZCart
+ * _ Four estimated components : Mean of XCart, YCart, ZCart and an isData mask
+ * _ Optionnal image : No
+ *
+ * \ingroup DiapOTBModule
+ */
+template <class TInputPixel, class TOutputPixel>
+class SARCartesianMeanFunctor : public SARPolygonsFunctor<TInputPixel, TOutputPixel>
+{
+public:
+  /** Standard class typedefs */
+  using Self         = SARCartesianMeanFunctor;
+  using Superclass   = SARPolygonsFunctor<TInputPixel, TOutputPixel>;
+  using Pointer      = itk::SmartPointer<Self>;
+  using ConstPointer = itk::SmartPointer<const Self>;
+
+  /** Runtime information */
+  itkTypeMacro(SARCartesianMeanFunctor, itk::Object);
+
+  /**
+   * Method SetRequiredComponentsToInd (override)
+   */
+  void SetRequiredComponentsToInd(std::map<std::string, int> const& mapForInd) final
+  {
+    std::map<std::string, int>::const_iterator it = mapForInd.begin();
+    while(it!=mapForInd.end())
+    {
+      if (it->first == "L")
+      {
+        m_indL = it->second;
+      }
+      else if (it->first == "C")
+      {
+        m_indC = it->second;
+      }
+      else if (it->first == "Y")
+      {
+        m_indY = it->second;
+      }
+      else if (it->first == "Z")
+      {
+        m_indZ = it->second;
+      }
+      else if (it->first == "XCart")
+      {
+        m_indXCart = it->second;
+      }
+      else if (it->first == "YCart")
+      {
+        m_indYCart = it->second;
+      }
+      else if (it->first == "ZCart")
+      {
+        m_indZCart = it->second;
+      }
+      ++it;
+    }
+  }
+
+  /**
+   * Method initialize (override)
+   */
+  void initalize(otb::Span<TOutputPixel> outValue, int threadId=1) final
+  {
+    int const nbComponentsMax = this->GetNumberOfEstimatedComponents();
+    int const nbElt = outValue.size();
+
+    m_CountPolygons[threadId].resize(nbElt);
+
+    std::fill(std::begin(m_CountPolygons[threadId]), std::end(m_CountPolygons[threadId]), 0);
+
+    // Reserve Components into outValue and Initalize
+    for(int ind = 0; ind < nbElt; ind++)
+    {
+      if (outValue[ind].GetSize() != nbComponentsMax)
+      {
+        outValue[ind].Reserve(nbComponentsMax);
+      }
+
+      for (int iC = 0; iC < nbComponentsMax; iC++)
+      {
+        outValue[ind][iC] = 0;
+      }
+    }
+  }
+
+  /**
+   * Method estimate (override)
+   */
+  void estimate(const int ind_LineSAR,
+      TInputPixel const& CLZY_InSideUp,
+      TInputPixel const& CLZY_InSideDown,
+      TInputPixel const& CLZY_OutSideUp,
+      TInputPixel const& CLZY_OutSideDown,
+      int firstCol_into_outputRegion,
+      otb::Span<TOutputPixel> outValue,
+      bool isFirstMaille, double & tgElevationMaxForCurrentLine,int threadId=1) final
+  {
+    // ML factor in range (to get SAR geometry (= Projeted DEM Pixel/Bands Geometry))
+    firstCol_into_outputRegion = m_originC + (firstCol_into_outputRegion - m_originC)*m_MLran;
+    int nbCol_into_outputRegion = outValue.size() * m_MLran;
+
+    if (static_cast<unsigned int>(nbCol_into_outputRegion+firstCol_into_outputRegion) >
+        (this->GetNbColSAR() + m_originC))
+    {
+      nbCol_into_outputRegion = this->GetNbColSAR() - firstCol_into_outputRegion + m_originC;
+    }
+
+    assert(threadId < m_CountPolygons.size());
+    auto & polygons_count = m_CountPolygons[threadId];
+
+    // Elevation angle (for shadow)
+    double tgElevationMax = 0.;
+
+
+    // Define the coef a and b
+    double a1 =  CLZY_InSideUp[m_indL] - ind_LineSAR;
+    a1 /= CLZY_InSideUp[m_indL] - CLZY_InSideDown[m_indL];
+
+    double const b1 = 1-a1;
+
+    double const c1 = b1*CLZY_InSideUp[m_indC] + a1*CLZY_InSideDown[m_indC];
+
+    double a2 = CLZY_OutSideUp[m_indL] - ind_LineSAR;
+    a2 /= (CLZY_OutSideUp[m_indL] - CLZY_OutSideDown[m_indL]);
+
+    double const b2 = 1-a2;
+
+    double const c2 = b2*CLZY_OutSideUp[m_indC] + a2*CLZY_OutSideDown[m_indC];
+
+    // Caculate zE, yE, zS and yS
+    // Select input and output side for the current maille with c1 and c2 comparison
+    double zE = b2*CLZY_OutSideUp[m_indZ] + a2*CLZY_OutSideDown[m_indZ];
+    double yE = b2*CLZY_OutSideUp[m_indY] + a2*CLZY_OutSideDown[m_indY];
+
+    double xCartE = b2*CLZY_OutSideUp[m_indXCart] + a2*CLZY_OutSideDown[m_indXCart];
+    double yCartE = b2*CLZY_OutSideUp[m_indYCart] + a2*CLZY_OutSideDown[m_indYCart];
+    double zCartE = b2*CLZY_OutSideUp[m_indZCart] + a2*CLZY_OutSideDown[m_indZCart];
+
+    double zS = b1*CLZY_InSideUp[m_indZ] + a1*CLZY_InSideDown[m_indZ];
+    double yS = b1*CLZY_InSideUp[m_indY] + a1*CLZY_InSideDown[m_indY];
+
+    double xCartS = b1*CLZY_InSideUp[m_indXCart] + a1*CLZY_InSideDown[m_indXCart];
+    double yCartS = b1*CLZY_InSideUp[m_indYCart] + a1*CLZY_InSideDown[m_indYCart];
+    double zCartS = b1*CLZY_InSideUp[m_indZCart] + a1*CLZY_InSideDown[m_indZCart];
+
+    double cE = c2;
+    double cS = c1;
+
+    if (c2 >= c1)
+    {
+      std::swap(zE, zS);
+      std::swap(yE, yS);
+
+      std::swap(xCartE, xCartS);
+      std::swap(yCartE, yCartS);
+      std::swap(zCartE, zCartS);
+      std::swap(cE, cS);
+    }
+
+    double const dc = cS - cE;
+
+    // Colunm into // TODO: the current maille
+    int const col1 = int(std::min(cE, cS))+1;        // +1 ?
+    int const col2 = int(std::max(cE, cS));
+
+    // tgElevationMax per mailles (not sure)
+    if (isFirstMaille)
+    {
+      tgElevationMax = (zE > 0.0) ?  yE / zE : sign (1e20, yE);
+      if (tgElevationMax < 0.) tgElevationMax = 0.;
+    }
+    else
+    {
+      tgElevationMax = tgElevationMaxForCurrentLine;
+    }
+
+    // Loop on colunm
+    for (int k = std::max(col1, firstCol_into_outputRegion),
+        N = std::min(col2, nbCol_into_outputRegion+firstCol_into_outputRegion-1)
+        ; k <= N; ++k)
+    {
+      double const a = (cS - k) / dc;
+      double const b = 1.0 - a;
+      double const Z = a * zE + b * zS;
+      double const Y = a * yE + b * yS;
+      double const XCart = a*xCartE + b*xCartS;
+      double const YCart = a*yCartE + b*yCartS;
+      double const ZCart = a*zCartE + b*zCartS;
+
+      double const tgElevation = (Z > 0.0) ? Y / Z : sign (1e20, Y);
+
+      // Check into shadows?
+      bool const not_in_shadows = tgElevation > tgElevationMax;
+      int const validity_flag = not_in_shadows ? 1 : 2;
+      if (not_in_shadows)
+        tgElevationMax = tgElevation;
+
+      // Contribution for this maille
+      int ind_for_out = k - firstCol_into_outputRegion;
+
+      // Into ML Geo if m_MLRan > 1
+      ind_for_out /= m_MLran; // Integer division
+
+      // Accumulate X, Y and Z Cartesian
+      assert(ind_for_out < outValue.size());
+      outValue[ind_for_out][0] += XCart;
+      outValue[ind_for_out][1] += YCart;
+      outValue[ind_for_out][2] += ZCart;
+
+      // Set isData at 1
+      outValue[ind_for_out][3] = validity_flag;
+
+      // Increment the counter
+      assert(ind_for_out < polygons_count.size());
+      polygons_count[ind_for_out] += 1;
+    }
+
+    // Store the tgElevationMax (for the next mailles)
+    tgElevationMaxForCurrentLine = tgElevationMax;
+  }
+
+  /**
+   * Synthetize to apply on `XCart,` `YCart `and `ZCart` the counter of polygons.
+   */
+  void synthetize(otb::Span<TOutputPixel> outValue, int threadId=1) final
+  {
+    assert(threadId < m_CountPolygons.size());
+    auto const& polygons_count = m_CountPolygons[threadId];
+    for(int ind = 0, nbElt = outValue.size(); ind < nbElt; ind++)
+    {
+      assert(ind < polygons_count.size());
+      if (polygons_count[ind] !=0)
+      {
+        outValue[ind][0] /= polygons_count[ind];
+        outValue[ind][1] /= polygons_count[ind];
+        outValue[ind][2] /= polygons_count[ind];
+      }
+    }
+  }
+
+  // default constructor
+  SARCartesianMeanFunctor() = delete;
+
+  /** Constructor */ //
+  SARCartesianMeanFunctor(int nbThreads, unsigned int mlran = 1, unsigned int mlazi = 1,
+      unsigned int originC = 0)
+  : m_NbThreads(nbThreads)
+  , m_CountPolygons(m_NbThreads) // Allocate the buffer of polygons for each thread
+  , m_MLran( mlran )
+  , m_MLazi( mlazi )
+  , m_originC( originC )
+  {
+    // Global argument (from PolygonsFunctor)
+    this->SetNumberOfExpectedComponents(7); // expected components in input = 7
+    this->SetNumberOfEstimatedComponents(4); // estimated components in output = 4
+
+    if (m_MLran == 1 && m_MLazi == 1)
+    {
+      this->SetOutputGeometry("SAR");
+    }
+    else
+    {
+      this->SetOutputGeometry("ML");
+    }
+
+    // Set the vector of Required Components
+    std::vector<std::string> vecComponents{
+      "C", "L", "Z", "Y", "XCart", "YCart", "ZCart"};
+    this->SetRequiredComponents(move(vecComponents));
+
+    std::vector<std::string> estimatedComponents{"XCart", "YCart", "ZCart", "isData"};
+    this->SetEstimatedComponents(move(estimatedComponents));
+  }
+
+  using Superclass::GetOutputGeometry;
+
+  /**
+   * Method GetOutputGeometry (override)
+   */
+  void GetOutputGeometry(std::string & outputGeo, unsigned int & mlran, unsigned int & mlazi) final
+  {
+    if (m_MLran == 1 && m_MLazi == 1)
+    {
+      outputGeo = "SAR";
+    }
+    else
+    {
+      outputGeo = "ML";
+    }
+
+    mlran = m_MLran;
+    mlazi = m_MLazi;
+  }
+
+  // Redefinition to use construction with NumberOfThreads and ML factor (Functor used in Multi-Threading)
+  static Pointer New(int nbThreads, unsigned int mlran = 1, unsigned int mlazi = 1, unsigned int originC = 0)
+  {
+    Pointer smartPtr = ::itk::ObjectFactory< Self >::Create();
+    if ( smartPtr == nullptr )
+    {
+      smartPtr = new Self(nbThreads, mlran, mlazi, originC);
+    }
+    smartPtr->UnRegister();
+    return smartPtr;
+  }
+
+  /** Destructor */
+  ~SARCartesianMeanFunctor() final = default;
+
+protected :
+  int m_NbThreads;
+  std::vector<std::vector<int>> m_CountPolygons;
+
+  // Index for each components
+  int m_indL;
+  int m_indC;
+  int m_indZ;
+  int m_indY;
+  int m_indXCart;
+  int m_indYCart;
+  int m_indZCart;
+
+  // ML factors
+  unsigned int m_MLran;
+  unsigned int m_MLazi;
+
+  // Origin for range dimension
+  unsigned int m_originC;
+};
+
+}
+}
+
+#endif

include/otbSARDEMPolygonsAnalysisImageFilter2.h

+/*
+ * 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 otbSARDEMPolygonsAnalysisImageFilter2_h
+#define otbSARDEMPolygonsAnalysisImageFilter2_h
+
+#include "otbSpan.h"
+#include "otbImageKeywordlist.h"
+#include "otbGenericRSTransform.h"
+
+#include "itkImageToImageFilter.h"
+#include "itkSmartPointer.h"
+#include "itkPoint.h"
+
+#include "itkImageScanlineConstIterator.h"
+#include "itkImageScanlineIterator.h"
+#include "itkImageRegionIterator.h"
+
+#include <type_traits>
+
+namespace otb
+{
+/** \class SARDEMPolygonsAnalysisImageFilter
+ * \brief Analyses the projected DEM.
+ *
+ * This filter :
+ * _ builds polygons with projected DEM points.
+ * _ select DEM polygons thanks to intersection with each SAR line.
+ * _ scan all selected polygons to extract wanted information for each pixel into output geometry. The "wanted
+ * information" is defined by the TFunction. This functor has to be instanciated outside this filter and is set
+ * with SetFunctorPtr.
+ *
+ * This filter can be used to estimate the amplitude image of a SAR image. Two kinds of output geometries are
+ * available : SAR geometry (same as TImageSAR) or ML geometry (thanks to ML factors). A optional image defined
+ * as a vector can be estimated.
+ *
+ * The output geometry and the optional image are required by the functor itself. Outputs can be a simple images
+ * or vector images.
+ *
+ * \ingroup DiapOTBModule
+ */
+
+template <typename TImageIn,  typename TImageOut, typename TImageDEM, typename TImageSAR, typename TFunction>
+class ITK_EXPORT SARDEMPolygonsAnalysisImageFilter
+: public itk::ImageToImageFilter<TImageIn,TImageOut>
+{
+public:
+
+  // Standard class typedefs
+  using Self         = SARDEMPolygonsAnalysisImageFilter;
+  using Superclass   = itk::ImageToImageFilter<TImageIn,TImageOut>;
+  using Pointer      = itk::SmartPointer<Self>;
+  using ConstPointer = itk::SmartPointer<const Self>;
+
+  // Method for creation through object factory
+  itkNewMacro(Self);
+  // Run-time type information
+  itkTypeMacro(SARDEMPolygonsAnalysisImageFilter,ImageToImageFilter);
+
+  /** Typedef to image input type : otb::VectorImage */
+  using ImageInType                = TImageIn;
+  /** Typedef to describe the inout image pointer type. */
+  using ImageInPointer             = typename ImageInType::Pointer;
+  using ImageInConstPointer        = typename ImageInType::ConstPointer;
+  /** Typedef to describe the inout image region type. */
+  using ImageInRegionType          = typename ImageInType::RegionType;
+  /** Typedef to describe the type of pixel and point for inout image. */
+  using ImageInPixelType           = typename ImageInType::PixelType;
+  using ImageInPointType           = typename ImageInType::PointType;
+  using ImageInSubPixelType        = typename ImageInType::InternalPixelType;
+  static_assert(std::is_floating_point<ImageInSubPixelType>::value, "float or double!");
+  using FastPixelType              = ImageInSubPixelType const*;
+
+  /** Typedef to describe the image index, size types and spacing for inout image. */
+  using ImageInIndexType           = typename ImageInType::IndexType;
+  using ImageInIndexValueType      = typename ImageInType::IndexValueType;
+  using ImageInSizeType            = typename ImageInType::SizeType;
+  using ImageInSizeValueType       = typename ImageInType::SizeValueType;
+  using ImageInSpacingType         = typename ImageInType::SpacingType;
+  using ImageInSpacingValueType    = typename ImageInType::SpacingValueType;
+
+  /** Typedef to image output type : otb::Image or otb::VectorImage with SAR or ML geometry  */
+  using ImageOutType               = TImageOut;
+  /** Typedef to describe the output image pointer type. */
+  using ImageOutPointer            = typename ImageOutType::Pointer;
+  using ImageOutConstPointer       = typename ImageOutType::ConstPointer;
+  /** Typedef to describe the output image region type. */
+  using ImageOutRegionType         = typename ImageOutType::RegionType;
+  /** Typedef to describe the type of pixel and point for output image. */
+  using ImageOutPixelType          = typename ImageOutType::PixelType;
+  using ImageOutPointType          = typename ImageOutType::PointType;
+  /** Typedef to describe the image index, size types and spacing for output image. */
+  using ImageOutIndexType          = typename ImageOutType::IndexType;
+  using ImageOutIndexValueType     = typename ImageOutType::IndexValueType;
+  using ImageOutSizeType           = typename ImageOutType::SizeType;
+  using ImageOutSizeValueType      = typename ImageOutType::SizeValueType;
+  using ImageOutSpacingType        = typename ImageOutType::SpacingType;
+  using ImageOutSpacingValueType   = typename ImageOutType::SpacingValueType;
+
+  /** Typedef to optionnal image output type : otb::Image or otb::VectorImage. This output is a vector
+   * with dimensions 1 x nbLineOfMainOutput */
+  using ImageOptionnalType         = TImageOut;
+  /** Typedef to describe the optionnal output image pointer type. */
+  using ImageOptionnalPointer      = typename ImageOptionnalType::Pointer;
+  using ImageOptionnalConstPointer = typename ImageOptionnalType::ConstPointer;
+  using ImageOptionnalRegionType   = typename ImageOptionnalType::RegionType;
+  using ImageOptionnalIndexType    = typename ImageOptionnalType::IndexType;
+  using ImageOptionnalSizeType     = typename ImageOptionnalType::SizeType;
+
+
+  /** Typedef to DEM image type : otb::Image with DEM geometry (only metadata)  */
+  using ImageDEMType               = TImageDEM;
+  using ImageDEMPointer            = typename ImageDEMType::Pointer;
+
+  /** Typedef to SAR image type : otb::Image with SAR geometry (only metadata)  */
+  using ImageSARType               = TImageSAR;
+  using ImageSARPointer            = typename ImageSARType::Pointer;
+
+  // Define Point2DType and Point3DType
+  using Point2DType = itk::Point<double,2>;
+  using Point3DType = itk::Point<double,3>;
+
+  // Define RSTRansform For Localisation inverse (SAR to DEM)
+  using RSTransformType2D = otb::GenericRSTransform<double,2,2>;
+
+  using FunctionType      = TFunction;
+  using FunctionPointer   = typename FunctionType::Pointer;
+
+  // Typedef for iterators
+  using InputIterator     = itk::ImageScanlineConstIterator< ImageInType >;
+  using OutputIterator    = itk::ImageScanlineIterator< ImageOutType >;
+  using OutputOptIterator = itk::ImageRegionIterator< ImageOutType >;
+
+  // Setter
+  void SetSARImageKeyWorList(ImageKeywordlist const& sarImageKWL);
+  void SetSARImagePtr(ImageSARPointer sarPtr);
+  void SetDEMImagePtr(ImageDEMPointer demPtr);
+  void SetDEMInformation(double gain, int DEMDirC, int DEMDirL);
+  void SetFunctorPtr(FunctionPointer functor);
+  itkSetMacro(NoData, int);
+  // Getter
+  itkGetConstMacro(Gain, double);
+  itkGetMacro(NoData, int);
+
+  // Initialize margin and GenericRSTransform
+  void initializeMarginAndRSTransform();
+
+protected:
+  // Default Constructor
+  SARDEMPolygonsAnalysisImageFilter() = default;
+
+  // Destructor
+  ~SARDEMPolygonsAnalysisImageFilter() final = default;
+
+  // Print
+  void PrintSelf(std::ostream & os, itk::Indent indent) const final;
+
+  /**
+   * SARDEMPolygonsAnalysisImageFilter produces an image which are into SAR geometry or ML geometry.
+   * The differences between output and input images are the size of images and the dimensions. The input image
+   * are a otb::VectorImage and the output can be a classic image or a vector image.
+   * As such, SARDEMPolygonsAnalysisImageFilter needs to provide an implementation for
+   * GenerateOutputInformation() in order to inform the pipeline execution model.
+   */
+  void GenerateOutputInformation() final;
+
+  /**
+   * SARDEMPolygonsAnalysisImageFilter needs a input requested region that corresponds to the projection
+   * into our main output requested region.
+   * As such, SARQuadraticAveragingImageFilter needs to provide an implementation for
+   * GenerateInputRequestedRegion() in order to inform the pipeline execution model.
+   * \sa `ProcessObject::GenerateInputRequestedRegion()`
+   */
+  void GenerateInputRequestedRegion() final;
+
+  /**
+   * OutputRegionToInputRegion returns the input region. This input region corresponds to the outputRegion.
+   */
+  ImageInRegionType OutputRegionToInputRegion(const ImageOutRegionType& outputRegion) const;
+
+  void BeforeThreadedGenerateData() final;
+
+  /**
+   * SARDEMPolygonsAnalysisImageFilter can be implemented as a multithreaded filter.
+   * Therefore, this implementation provides a ThreadedGenerateData() routine
+   * which is called for each processing thread. The main output image data is
+   * allocated automatically by the superclass prior to calling
+   * ThreadedGenerateData().  ThreadedGenerateData can only write to the
+   * portion of the output image specified by the parameter
+   * "outputRegionForThread"
+   *
+   * \sa `ImageToImageFilter::ThreadedGenerateData()`,
+   *     `ImageToImageFilter::GenerateData()`
+   */
+  void ThreadedGenerateData(const ImageOutRegionType& outputRegionForThread, itk::ThreadIdType threadId) final;
+
+  /**
+   * Select for the current SAR Line (ind_LineSAR), the polygons that intersect the current SAR line.
+   * This selection is always made with the SAR Line (even in ML geometry for the main output).
+   */
+  bool PolygonAndSarLineIntersection(const int ind_LineSAR,
+      const ImageInSubPixelType ligUL, const ImageInSubPixelType colUL,
+      const ImageInSubPixelType ligUR, const ImageInSubPixelType colUR,
+      const ImageInSubPixelType ligLR, const ImageInSubPixelType colLR,
+      const ImageInSubPixelType ligLL, const ImageInSubPixelType colLL,
+      int & i_InSideUp, int & i_InSideDown,
+      int & i_OutSideUp, int & i_OutSideDown);
+
+  /**
+   * Scan all selected polygons and apply on it, the choosen functor.
+   */
+  void PolygonsScan(const int ind_LineSAR,
+      std::vector<std::vector<FastPixelType> const*> const& CLZY_InSideUp_Vec,
+      std::vector<std::vector<FastPixelType> const*> const& CLZY_InSideDown_Vec,
+      std::vector<std::vector<FastPixelType> const*> const& CLZY_OutSideUp_Vec,
+      std::vector<std::vector<FastPixelType> const*> const& CLZY_OutSideDown_Vec,
+      int firstCol_into_outputRegion,
+      otb::Span<ImageOutPixelType> outValueTab, int threadId);
+
+private:
+  SARDEMPolygonsAnalysisImageFilter(const Self&) = delete;
+  Self& operator=(const Self &) = delete;
+
+  // SAR Image (only metadata)
+  ImageSARPointer m_SarImagePtr;
+
+  // DEM Image (only metadata)
+  ImageDEMPointer m_DemImagePtr;
+
+  // SAR Image KeyWorldList
+  ImageKeywordlist m_SarImageKwl;
+
+  // Multiplying gain to obtain a mean radiometry at 100
+  double m_Gain = 100.;
+
+  // Directions to scan the DEM in function of SAR geometry (from near to far range)
+  int m_DEMdirL;
+  int m_DEMdirC;
+
+  // Margin for polygon selection
+  int m_Margin = 0;
+
+  // Value for NoDATA
+  int m_NoData = -32768;
+
+  // Function to apply on DEM Polygon
+  FunctionPointer m_FunctionOnPolygon;
+
+  // Output Geometry (SAR or ML)
+  std::string m_OutGeometry;
+  // ML factor (for ML geometry)
+  unsigned int m_MLRan = 1;
+  unsigned int m_MLAzi = 1;
+
+  // Sar dimensions
+  int m_NbLinesSAR;
+  int m_NbColSAR;
+
+  // RSTransform (for inverse localisation)
+  RSTransformType2D::Pointer m_RSTransform;
+
+  std::vector<bool> m_is_lineX_sorted;
+  bool              m_all_lines_are_sorted = false;
+};
+
+} // End namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSARDEMPolygonsAnalysisImageFilter2.hxx"
+#endif
+
+#endif

include/otbSARPolygonsFunctor2.h

+/*
+ * 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 otbSARPolygonsFunctor2_h
+#define otbSARPolygonsFunctor2_h
+
+#include "otbSpan.h"
+#include "itkObject.h"
+
+#include <map>
+#include <string>
+#include <cmath>
+#include <algorithm>
+
+namespace
+{ // Anonymous namespace
+template <typename F>
+inline
+F sign(F a, F b)
+{
+  return b >= 0
+    ? std::fabs(a)
+    :-std::fabs(a);
+}
+} // Anonymous namespace
+
+namespace otb
+{
+namespace Function
+{
+/** \class SARPolygonsFunctor
+ * \brief Base class for functor used on polygons into SARDEMPolygonsAnalysis
+ *
+ * This base class contains four main functions:
+ * - initialize : to initialize the outValue array (one for each output line). Initialize a simple array by default.
+ * - estimate : to calculate the outValue array. Has to be implemented.
+ * - synthetize : to synthetize the estimation made into estimate function.
+ *
+ * Useful to calculate mean, for example. Empty By default.
+ * - estimateOptionnalImage : to estimate the optionnal image. Empty By default.
+ *
+ * Several setter/getter are available to specify each functor.
+ * For example each functor can choose :
+ *
+ * - The main output geometry (SAR or ML)
+ * - The required components into projeted DEM
+ * - The estimated components
+ * - The optionnal image estimation
+ *
+ * \sa PolygonsFunctor
+ *
+ * \ingroup DiapOTBModule
+ */
+template <class TInputPixel, class TOutputPixel>
+class SARPolygonsFunctor : public itk::Object
+{
+public:
+  /** Standard class typedefs */
+  using Self         = SARPolygonsFunctor;
+  using Superclass   = itk::Object;
+  using Pointer      = itk::SmartPointer<Self>;
+  using ConstPointer = itk::SmartPointer<const Self>;
+
+  /** Runtime information */
+  itkTypeMacro(SARPolygonsFunctor, itk::Object);
+
+  // Initialize method
+  virtual void initalize(otb::Span<TOutputPixel> outValue, int /*threadId=1*/)
+  {
+    // Re Initialize outValue
+    std::fill(std::begin(outValue), std::end(outValue), TOutputPixel{});
+  }
+
+  // Estimate method
+  virtual void estimate(const int ind_LineSAR,
+      TInputPixel const& CLZY_InSideUp,
+      TInputPixel const& CLZY_InSideDown,
+      TInputPixel const& CLZY_OutSideUp,
+      TInputPixel const& CLZY_OutSideDown,
+      int firstCol_into_outputRegion,
+      otb::Span<TOutputPixel> outValue,
+      bool isFirstMaille, double & tgElevationMaxForCurrentLine, int threadId=1) = 0;
+
+  // Synthetize method
+  virtual void synthetize(otb::Span<TOutputPixel> , int /*threadId=1*/) = 0;
+
+  // estimateOptionnalImage method
+  virtual void estimateOptionnalImage(TOutputPixel * /*outValue*/, int /*threadId=1*/) {}
+
+  // Getter/Setter
+  itkSetMacro(NumberOfExpectedComponents, unsigned int);
+  itkSetMacro(NumberOfEstimatedComponents, unsigned int);
+  itkSetMacro(OutputGeometry, std::string);
+  itkSetMacro(WithOptionnalOutput, bool);
+
+  void SetEstimatedComponents(std::vector<std::string> vecEstimated)
+  {
+    m_EstimatedComponents = move(vecEstimated);
+  }
+
+  void SetRequiredComponents(std::vector<std::string> vecRequired)
+  {
+    m_RequiredComponents = move(vecRequired);
+  }
+
+  void SetSARDimensions(int nbColSAR, int nbLinesSAR)
+  {
+    m_NbColSAR = nbColSAR;
+    m_NbLinesSAR = nbLinesSAR;
+  }
+
+  virtual void SetRequiredComponentsToInd(std::map<std::string, int> const& mapForInd)
+  {
+    std::map<std::string, int>::const_iterator it = mapForInd.begin();
+    while(it!=mapForInd.end())
+    {
+      m_RequiredComponentsToInd[it->first] = it->second;
+      ++it;
+    }
+  }
+
+  itkGetMacro(NumberOfExpectedComponents, unsigned int);
+  itkGetMacro(NumberOfEstimatedComponents, unsigned int);
+  itkGetMacro(EstimatedComponents, std::vector<std::string>);
+  itkGetMacro(RequiredComponents, std::vector<std::string>);
+  itkGetMacro(OutputGeometry, std::string);
+  itkGetMacro(NbColSAR, int);
+  itkGetMacro(NbLinesSAR, int);
+  itkGetMacro(WithOptionnalOutput, bool);
+  virtual void GetOutputGeometry(std::string & outputGeo, unsigned int & mlran, unsigned int & mlazi)
+  {
+    outputGeo = m_OutputGeometry;
+    mlran = 1;
+    mlazi = 1;
+  }
+  std::map<std::string, int> & GetRequiredComponentsToInd() const
+  {
+    return m_RequiredComponentsToInd;
+  }
+
+  /** Default constructor */
+  SARPolygonsFunctor() = default;
+
+  /** Destructor */
+  ~SARPolygonsFunctor() override = default;
+
+private :
+
+  // Number of expected Components for input
+  unsigned m_NumberOfExpectedComponents = 1;
+  // Number of estimated Components for output
+  unsigned m_NumberOfEstimatedComponents = 1;
+  // Geometry of output image into PolygonsAnalysis Filter (SAR or Polygons)
+  std::string m_OutputGeometry = "SAR";
+  // Vector of required Components
+  std::vector<std::string> m_RequiredComponents;
+  // Map of required Components
+  std::map<std::string, int> m_RequiredComponentsToInd;
+  // Vector of estimated Components
+  std::vector<std::string> m_EstimatedComponents;
+  // SAR Dimensions
+  int m_NbColSAR;
+  int m_NbLinesSAR;
+  // Boolean for optionnal output
+  bool m_WithOptionnalOutput = false; // Optionnal output must be a complement of the main output (small image)
+};
+
+}
+}
+
+#endif

include/otbVLVPointIterator.h

+/*
+ * 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 otbVLVPointVLVPointIterator_h
+#define otbVLVPointVLVPointIterator_h
+
+#include "itkVariableLengthVector.h"
+#include <iterator>
+#include <type_traits>
+#include <cassert>
+
+namespace otb
+{
+
+/**
+ * C++ standard compliant iterator for contiguous arrays of
+ * `itk::VariableLengthVector`.
+ * This minimalist iterator permits the building block to define iterable
+ * ranges over set of `itk::VariableLengthVector` that can then be used in C++
+ * standards algorithms.
+ *
+ * \internal The important part is that `+1` needs to take the number of
+ * components/bands into account.
+ * \tparam T  sub-pixel data type
+ */
+template <typename T>
+struct VLVPointIterator
+{
+  using iterator_category = std::random_access_iterator_tag;
+  using value_type        = T;
+  using difference_type   = std::ptrdiff_t;
+  using pointer           = T*;
+  using reference         = itk::VariableLengthVector<typename std::remove_cv<T>::type>;
+
+  explicit VLVPointIterator(T* p, unsigned nb_bands) noexcept
+    : m_p(p), m_nb_bands(nb_bands)
+    {}
+
+  VLVPointIterator& operator+=(difference_type o) noexcept {
+    m_p += o * m_nb_bands;
+    return *this;
+  }
+  friend VLVPointIterator operator+(VLVPointIterator it, difference_type o) noexcept {
+    it += o;
+    return it;
+  }
+  friend VLVPointIterator operator+(difference_type o, VLVPointIterator it) noexcept {
+    it += o;
+    return it;
+  }
+
+  VLVPointIterator& operator-=(difference_type o) noexcept {
+    m_p -= o * m_nb_bands;
+    return *this;
+  }
+  friend VLVPointIterator operator-(VLVPointIterator it, difference_type o) noexcept {
+    it -= o;
+    return it;
+  }
+
+  friend difference_type operator-(VLVPointIterator const& lhs, VLVPointIterator const& rhs) noexcept {
+    assert(lhs.m_nb_bands == rhs.m_nb_bands);
+    return (lhs.m_p - rhs.m_p) / lhs.m_nb_bands;
+  }
+
+  VLVPointIterator& operator++() noexcept {
+    m_p += m_nb_bands;
+    return *this;
+  }
+  VLVPointIterator operator++(int) noexcept {
+    auto tmp = *this;
+    m_p += m_nb_bands;
+    return tmp;
+  }
+
+  VLVPointIterator& operator--() noexcept {
+    m_p -= m_nb_bands;
+    return *this;
+  }
+  VLVPointIterator operator--(int) noexcept {
+    auto tmp = *this;
+    m_p -= m_nb_bands;
+    return tmp;
+  }
+
+#if 0
+  T& operator[](unsigned p) noexcept {
+    assert(p < m_nb_bands);
+    return m_p[p];
+  }
+  T const& operator[](unsigned p) const noexcept {
+    assert(p < m_nb_bands);
+    return m_p[p];
+  }
+#endif
+
+  reference operator*() const noexcept {
+    return itk::VariableLengthVector<typename std::remove_cv<T>::type>(m_p, m_nb_bands);
+  }
+
+  decltype(auto) data() const noexcept { return m_p; }
+
+  friend bool operator==(VLVPointIterator const& rhs, VLVPointIterator const& lhs) noexcept
+  { return rhs.m_p == lhs.m_p; }
+  friend bool operator!=(VLVPointIterator const& rhs, VLVPointIterator const& lhs) noexcept
+  { return rhs.m_p != lhs.m_p; }
+  friend bool operator<=(VLVPointIterator const& rhs, VLVPointIterator const& lhs) noexcept
+  { return rhs.m_p <= lhs.m_p; }
+  friend bool operator<(VLVPointIterator const& rhs, VLVPointIterator const& lhs) noexcept
+  { return rhs.m_p < lhs.m_p; }
+  friend bool operator>(VLVPointIterator const& rhs, VLVPointIterator const& lhs) noexcept
+  { return rhs.m_p > lhs.m_p; }
+  friend bool operator>=(VLVPointIterator const& rhs, VLVPointIterator const& lhs) noexcept
+  { return rhs.m_p >= lhs.m_p; }
+
+  // itk::VariableLengthVector<typename std::remove_cv<T>::type> operator*() noexcept {
+    // return itk::VariableLengthVector<T>(m_p, m_nb_bands);
+  // }
+
+private:
+  T*       m_p;
+  unsigned m_nb_bands;
+};
+
+} // otb namespace
+
+#endif  // otbVLVPointVLVPointIterator_h