otbNormalCompute.h

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

#ifndef otbNormalCompute_h
#define otbNormalCompute_h

#include <vnl/vnl_cross.h>
#include <tuple>
#include <cmath>

namespace otb
{

/**
 * Helper class to define strong/opaque types.
 * The idea is to prevent easy mismatches whe passing multiples parameters of
 * the sam type.
 * \tparam T  Encapsulated type
 * \tparam Tag  Tag type used to distinguish different strong types.
 */
template <typename T, typename Tag>
struct StrongType
{
  StrongType() = default;
  explicit StrongType(T v) : value{v} {}
  explicit operator T const&() const noexcept { return value; }
private:
  T value;
};

namespace Tags
{
  struct Above{};
  struct Below{};
  struct Left {};
  struct Right{};
} // namespace Tags

template <typename T> using Above = StrongType<T, Tags::Above>;
template <typename T> using Below = StrongType<T, Tags::Below>;
template <typename T> using Left  = StrongType<T, Tags::Left>;
template <typename T> using Right = StrongType<T, Tags::Right>;

/**
 * Returns the lambda that'll be used to produce the normal vector at the
 * center of the neighbourhood.
 *
 * \tparam RealType meant to be `float` or `double`
 * \param[in] x_spacing spacing in X
 * \param[in] y_spacing spacing in Y
 * \pre spacings shall be expressed in the same unit as the heights.
 *
 * This flavour supposes a constant spacing between height _pixels_.
 *
 * If we note:
 * - \f$o\f$ the center
 * - \f$a\f$ the pixel above
 * - \f$b\f$ the pixel on the right
 * - \f$c\f$ the pixel below
 * - \f$d\f$ the pixel on the left
 * - \f$X\f$ the x-axis spacing
 * - \f$Y\f$ the y-axis spacing
 * - \f$H_{point}\f$ the height at the given _point_
 *
 * We can compute the normal vector of
 * - \f$\overrightarrow{ob} \wedge \overrightarrow{oa} = \begin{pmatrix}-Y.(H_b - H_o) \\ -X.(H_a - H_o) \\ X.Y\end{pmatrix}\f$
 * - \f$\overrightarrow{oc} \wedge \overrightarrow{ob} = \begin{pmatrix}-Y.(H_b - H_o) \\ +X.(H_c - H_o) \\ X.Y\end{pmatrix}\f$
 * - \f$\overrightarrow{od} \wedge \overrightarrow{oc} = \begin{pmatrix}+Y.(H_d - H_o) \\ +X.(H_c - H_o) \\ X.Y\end{pmatrix}\f$
 * - \f$\overrightarrow{oa} \wedge \overrightarrow{od} = \begin{pmatrix}+Y.(H_d - H_o) \\ -X.(H_a - H_o) \\ X.Y\end{pmatrix}\f$
 *
 * Their sum is \f$\begin{pmatrix}2.Y(H_d - H_b) \\ 2.X.(H_c-H_a) \\ 4.X.Y\end{pmatrix}\f$. The
 * lamdba used will return the normalized value.
 */
template <typename RealType>
inline
auto SimpleNormalMethodForRegularGrid(RealType x_spacing, RealType y_spacing)
{
  auto const z = std::abs(RealType{2} * x_spacing * y_spacing);
  auto const z2 = z * z;
  auto lambda = [=](
      Above<RealType> above, Left <RealType> left,
      Right<RealType> right, Below<RealType> below)
  {
    auto const x = y_spacing * (RealType(left)  - RealType(right));
    auto const y = x_spacing * (RealType(below) - RealType(above));
    auto const norm = 1.0f / std::sqrt(x*x + y*y + z2);

    return std::make_tuple(x*norm, y*norm, z*norm);

  };
  return lambda;
}

/**
 * Returns the lambda that'll be used to produce the normal vector at the
 * center of the neighbourhood.
 *
 * \tparam PointType Models {X, Y, Z} points
 *
 * This flavour expects point coordinates.
 *
 * If we note:
 * - \f$o\f$ the center
 * - \f$a\f$ the pixel above
 * - \f$b\f$ the pixel on the right
 * - \f$c\f$ the pixel below
 * - \f$d\f$ the pixel on the left
 *
 * We can compute the sum of all the normal vector to \f$<overrightarrow{ob}, \overrightarrow{oa}>\f$...
 * - \f$\overrightarrow{ob} \wedge \overrightarrow{oa}\f$
 * - \f$\overrightarrow{oc} \wedge \overrightarrow{ob}\f$
 * - \f$\overrightarrow{od} \wedge \overrightarrow{oc}\f$
 * - \f$\overrightarrow{oa} \wedge \overrightarrow{od}\f$
 *
 * Which, given cross product properties, is equal to
 * \f$\overrightarrow{ca} \wedge \overrightarrow{bd}\f$
 * The lamdba used will return the normalized value, or a triplet of `nodata`
 * if any of the two vertical or horizontal vector is null.
 */
template <typename PointType, typename FloatType>
inline
auto SimpleNormalMethodForCartesianPoints(FloatType nodata)
{
  auto lambda = [=](
      Above<PointType> above, Left <PointType> left,
      Right<PointType> right, Below<PointType> below)
  {
    auto const CA = PointType(above) - PointType(below);
    auto const BD = PointType(left)  - PointType(right);
    auto const n  = vnl_cross_3d(CA, BD);
    auto const norm = n.two_norm();
    if (norm == 0)
    {
      return PointType{nodata, nodata, nodata};
    }
    // assert(norm > 0);
    auto const inv_norm = static_cast<decltype(CA[0])>(1.0) / norm;

    return n * inv_norm;
  };
  return lambda;
}
} // otb namespace

#endif  // otbNormalCompute_h