otbExtractNormalVector.cxx

/*
 * Copyright (C) 2005-2020 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 "otbNormalCompute.h"
#include "otbRepack.h"
#include "otbPositionHelpers.h"

#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbFunctorImageFilter.h"
#include "otbDEMHandler.h"
#include "otbSarSensorModelAdapter.h"
// Elevation handler
#include "otbWrapperElevationParametersHandler.h"

#include "itkFixedArray.h"
#include <array>
#include <cmath>

namespace // anonymous
{

/**
 * Returns the lambda that'll be used to produce the normal vector at the
 * center of the neighbourhood.
 *
 * \note As a simplification, this method ignores Geiod correction on
 * neighbourhood pixels. They are all supposed identical.
 * \todo Handle cases with nodata => assume the data from the center; 0 if the
 * center is nodata
 */
template <typename OutputPixelType, typename InputImageType>
inline
auto AppSimpleNormalMethod(
    double            x_spacing,
    double            y_spacing,
    InputImageType  & image)
{
  using InputPixelType           = typename InputImageType::PixelType;
  using RealType                 = typename itk::PixelTraits<InputPixelType>::ValueType;
  using V3                       = vnl_vector_fixed<RealType, 3>;
  using NeighborhoodIteratorType = itk::ConstNeighborhoodIterator<InputImageType>;

  // x_spacing > 0 and y_spacing > 0 ==>
  //    6  7  8
  //    3  4  5
  //    0  1  2
  // x_spacing > 0 and y_spacing < 0 ==>
  //    0  1  2
  //    3  4  5
  //    6  7  8
  // x_spacing < 0 and y_spacing > 0 ==>
  //    8  7  6
  //    5  4  3
  //    2  1  0
  auto const above  = [y_spacing](NeighborhoodIteratorType const& n) { return otb::Above<V3>(otb::repack<V3>(n.GetPixel(y_spacing > 0 ? 7 : 1)));};
  auto const below  = [y_spacing](NeighborhoodIteratorType const& n) { return otb::Below<V3>(otb::repack<V3>(n.GetPixel(y_spacing > 0 ? 1 : 7)));};
  auto const left   = [x_spacing](NeighborhoodIteratorType const& n) { return otb::Left <V3>(otb::repack<V3>(n.GetPixel(x_spacing > 0 ? 3 : 5)));};
  auto const right  = [x_spacing](NeighborhoodIteratorType const& n) { return otb::Right<V3>(otb::repack<V3>(n.GetPixel(x_spacing > 0 ? 5 : 3)));};
  auto const center = [](NeighborhoodIteratorType const& n) { return n.GetPixel(4);};

  auto l0 = otb::SimpleNormalMethodForCartesianPoints<vnl_vector_fixed<double, 3>>();

  auto lambda = [=](NeighborhoodIteratorType const& in)
  {
    // Compute n-> at current point
    assert(in.Size() == 9);

    auto valid = [](auto v) /*-> decltype(v)*/ {
      // TODO: remove if we can't have nodata
      return v;
    };

    V3 const n = l0(above(valid(in)), valid(left(in)), valid(right(in)), valid(below(in)));
    static int i = 0;
    if (++i % 1000000 == 0)
    {
      std::cout << "-----@" << in.GetIndex() << " from " <<
        V3(above(in)) << V3(left(in)) << V3(right(in)) << V3(below(in))
        << " --> " << n
        <<std::endl;
    }

    auto const res = otb::repack<OutputPixelType>(n);
    return res;
  };
  return lambda;
}

} // anonymous namespace.

namespace otb
{

namespace Wrapper
{

/**
 * This application takes cares of building the normal vector from a DEM image.
 *
 * \author Luc Hermitte
 * \copyright CS Group
 * \ingroup App???
 * \ingroup SARCalibrationExtended
 */
class ExtractNormalVector : public Application
{
public:
  /** Standard class typedefs. */
  using Self             = ExtractNormalVector;
  using InputPixelType   = itk::FixedArray<double, 3>;
  using InputImageType   = otb::Image<InputPixelType, 2>;
  using OutputPixelType  = itk::VariableLengthVector<double>;

  // using InputImageType   = otb::Image<OutputPixelType, 2>;
  // using WriterType       = otb::ImageFileWriter<InputImageType>;

  /** Standard macro */
  itkNewMacro(Self);

  itkTypeMacro(ExtractNormalVector, unused);

private:
  void DoInit() override
  {
    SetName("ExtractNormalVector");
    SetDescription("This application produces normal vectors from a DEM image");

    // Documentation
    SetDocLongDescription("This application produces normal vectors from a DEM image" );


    SetDocLimitations("None");
    SetDocAuthors("Luc Hermitte (CS Group)");

    AddDocTag(Tags::Calibration);

    AddParameter(ParameterType_InputImage, "xyz",  "filename for the file containing XYZ cartesian coordinates aligned of their file geometry");
    SetParameterDescription("xyz", "filename for the file containing XYZ cartesian coordinates aligned of their file geometry");

    AddParameter(ParameterType_OutputImage, "out",  "Image of normal vectors");
    SetParameterDescription("out",  "Image of normal vectors");

    // AddParameter(ParameterType_Float, "nodata", "NoData value");
    // SetParameterDescription("nodata", "DEM empty cells are filled with this value (optional -32768 by default)");
    // SetDefaultParameterFloat("nodata", -32768);
    // MandatoryOff("nodata");

    AddRAMParameter();

    SetDocExampleParameterValue("xyz", "xyz.tif");
    SetDocExampleParameterValue("out", "normal.tif");
    SetOfficialDocLink();
  }

  void DoUpdateParameters() override
  {
  }

  void DoExecute() override
  {
    auto* inputImage = GetParameterImage<InputImageType>("xyz");
    assert(inputImage);
    // auto* inputImage = GetParameterFloatImage("dem");
    // auto* inputImage = GetParameterDoubleImage("dem");
    // auto* inputImage = GetParameterImage("dem");

    // auto const nodata = GetParameterFloat("nodata");
    // std::cout << "nodata      => " << nodata << std::endl;

    auto const spacing = inputImage->GetSignedSpacing();
    double const x_spacing = spacing[0];
    double const y_spacing = spacing[1];
    auto const lambda = AppSimpleNormalMethod<OutputPixelType>(x_spacing, y_spacing, *inputImage);

    auto filter = otb::NewFunctorFilter(lambda, 6, {{1,1}});
    filter->SetInput(inputImage);

    SetParameterOutputImage("out", filter->GetOutput());
    // SetParameterOutputImagePixelType("out",ImagePixelType_float);

    // Call register pipeline to allow streaming and garbage collection
    RegisterPipeline();
  }
};

} //end namespace Wrapper
} //end namespace otb

OTB_APPLICATION_EXPORT(otb::Wrapper::ExtractNormalVector)