Skip to content
Snippets Groups Projects
otbExtractNormalVector.cxx 8.79 KiB
Newer Older
Luc Hermitte's avatar
Luc Hermitte committed
 * Copyright (C) 2005-2021 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 "otbBandIds.h"
#include "otbRepack.h"
#include "otbPositionHelpers.h"

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

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

namespace // anonymous
{

using FloatType = float;

/**
 * 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(
    FloatType         x_spacing,
    FloatType         y_spacing,
    FloatType         nodata,
    InputImageType  & image)
  using InputPixelType           = typename InputImageType::PixelType;
  using RealType                 = typename itk::PixelTraits<InputPixelType>::ValueType;
  using InputV3                  = vnl_vector_fixed<RealType, 3>;
  using NeighborhoodIteratorType = itk::ConstNeighborhoodIterator<InputImageType>;

  auto const kwl = image.GetImageKeywordlist();
  auto const xband = value_or_unless(kwl, otb::bands::k_id_X, "extracting number of X band", 0);
  auto const yband = value_or_unless(kwl, otb::bands::k_id_Y, "extracting number of Y band", 1);
  auto const zband = value_or_unless(kwl, otb::bands::k_id_Z, "extracting number of Z band", 2);
  otbLogMacro(Info, <<"Using XYZ cartesian coordinates in bands " << xband << ", " << yband << ", " << zband);

  auto const XYZ = [xband, yband, zband](auto const& pixel) {
    return otb::repack<InputV3>(pixel[xband], pixel[yband], pixel[zband]);
  // 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_idx  = y_spacing > 0 ? 7 : 1;
  auto const below_idx  = y_spacing > 0 ? 1 : 7;
  auto const left_idx   = x_spacing > 0 ? 3 : 5;
  auto const right_idx  = x_spacing > 0 ? 5 : 3;
  auto const center_idx = 4;
  // TODO: Can we be sure this is enough to guarantee dot(res, center) > 0?
  auto const above  = [above_idx,  XYZ](NeighborhoodIteratorType const& n) { return otb::Above<InputV3>(XYZ(n.GetPixel(above_idx)));};
  auto const below  = [below_idx,  XYZ](NeighborhoodIteratorType const& n) { return otb::Below<InputV3>(XYZ(n.GetPixel(below_idx)));};
  auto const left   = [left_idx,   XYZ](NeighborhoodIteratorType const& n) { return otb::Left <InputV3>(XYZ(n.GetPixel(left_idx)));};
  auto const right  = [right_idx,  XYZ](NeighborhoodIteratorType const& n) { return otb::Right<InputV3>(XYZ(n.GetPixel(right_idx)));};
  auto const center = [center_idx, XYZ](NeighborhoodIteratorType const& n) { return XYZ(n.GetPixel(center_idx));};

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

  auto const nodata_res = otb::repack<OutputPixelType>(nodata, nodata, nodata);
  auto lambda = [=](NeighborhoodIteratorType const& in)
  {
    // Compute n-> at current point
    assert(in.Size() == 9);

    auto const is_invalid = [nodata](auto v)
    {
      InputV3 const& w = v; // Remove the StrongType<> layer
      return
        // SARCartesianMeanEstimation will return {0, 0, 0}
        w == InputV3{0., 0., 0.}
        // ossim may produce nan
        || std::isnan(w[0]) || std::isnan(w[1]) || std::isnan(w[2])
#if 0
        // other tools may produce nodata;
        // but we only support with SARDEMProjection => #if 0
        || w[0] == nodata || w[1] == nodata || w[2] == nodata
    auto const a = above(in);
    auto const l = left(in);
    auto const r = right(in);
    auto const b = below(in);
    auto const c = center(in);
    // TODO: Try to produce something if any left/right/above/below pixels is
    // valid => do a mean on the pixels of a side
    if (is_invalid(a) || is_invalid(l) || is_invalid(r) || is_invalid(b) ||
        is_invalid(c))
    {
      return nodata_res;
    }

    InputV3 n = l0(a, l, r, b);
    if (InputV3{0., 0., 0.} == n) // only way to be invalid => {0, 0, 0}
    // if (is_invalid(n))
    assert(! std::isnan(dot_product(n, c)));
    if (dot_product(n, c) < 0)
    {
      // TODO: we could expect this is the case for every single vector
      // computed for this image => factorize
      n = -n;
    }
    assert(dot_product(n, c) >= 0);
    if (++i % 10000000 == 0)
      std::cout << "-----@" << in.GetIndex() << " from " <<
        InputV3(above(in)) << " ; " << InputV3(left(in)) <<  " ; " << InputV3(right(in)) <<  " ; " << InputV3(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 ApplicationWithNoData
{
public:
  /** Standard class typedefs. */
  using Self             = ExtractNormalVector;
  // using InputPixelType   = itk::FixedArray<FloatType, 3>;
  /** Pixel Type of input image: VLV.
   * Some XYZ producers may provide mode than 3 bands... */
  using InputPixelType   = itk::VariableLengthVector<FloatType>;
  // using InputImageType   = otb::Image<InputPixelType, 2>;
  using InputImageType   = otb::VectorImage<FloatType, 2>;
  using OutputPixelType  = itk::VariableLengthVector<FloatType>;
  // 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");

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

  void DoUpdateParameters() override
  {
  }

  void DoExecute() override
  {
    auto* inputImage = GetParameterImage<InputImageType>("xyz");
    assert(inputImage);
    auto const nodata = GetParameterFloat("nodata");
    // std::cout << "nodata      => " << nodata << std::endl;
    AddNodataInMetadata(nodata);

    auto const spacing = inputImage->GetSignedSpacing();
    FloatType const x_spacing = spacing[0];
    FloatType const y_spacing = spacing[1];
    auto const lambda = AppSimpleNormalMethod<OutputPixelType>(x_spacing, y_spacing, nodata, *inputImage);
    auto filter = otb::NewFunctorFilter(lambda, 3, {{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)