Skip to content
Snippets Groups Projects

Draft: Resolve "Lot of variations in some places where topography varies a lot"

Merged Draft: Resolve "Lot of variations in some places where topography varies a lot"
3-lot-of-variations-in-some-places-where-topography-varies-a-lot into master
Merged Luc Hermitte requested to merge 3-lot-of-variations-in-some-places-where-topography-varies-a-lot into master
Viewing commit 1a6f5061
Next
Show latest version
1 file
+ 51
− 0
Compare changes
  • Side-by-side
  • Inline
app/otbExtractNormalVector.cxx
+ 214
− 52
@@ -35,44 +35,169 @@
#include "itkFixedArray.h"
#include <array>
#include <cmath>
#include <tuple>
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)
/** Type traits to factorize all type definitions. */
template <typename InputImageType>
struct input_types
{
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);
/** Helper functur to extract an {X, Y, Z} pixel from band information. */
template <typename InputImageType>
struct XYZ_t
{
using InputV3 = typename input_types<InputImageType>::InputV3;
auto const XYZ = [xband, yband, zband](auto const& pixel) {
explicit XYZ_t(InputImageType & image)
{
auto const kwl = image.GetImageKeywordlist();
xband = value_or_unless(kwl, otb::bands::k_id_X, "extracting number of X band", 0);
yband = value_or_unless(kwl, otb::bands::k_id_Y, "extracting number of Y band", 1);
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);
}
template <typename PixelType>
auto operator()(PixelType const& pixel) const {
return otb::repack<InputV3>(pixel[xband], pixel[yband], pixel[zband]);
};
}
private:
int xband;
int yband;
int zband;
};
template <typename InputImageType, bool does_compute_mean_on_sides>
struct side_extractors_t;
/** Specialization that returns a single pixel value for each side. */
template <typename InputImageType>
struct side_extractors_t<InputImageType, false>
{
side_extractors_t(
FloatType x_spacing,
FloatType y_spacing,
FloatType /*nodata*/,
InputImageType & image)
: XYZ(image)
, above_idx{y_spacing > 0 ? 7 : 1}
, below_idx{y_spacing > 0 ? 1 : 7}
, left_idx {x_spacing > 0 ? 3 : 5}
, right_idx{x_spacing > 0 ? 5 : 3}
{}
using NeighborhoodIteratorType = typename input_types<InputImageType>::NeighborhoodIteratorType;
using InputV3 = typename input_types<InputImageType>::InputV3;
auto above(NeighborhoodIteratorType const& n) const { return otb::Above<InputV3>(XYZ(n.GetPixel(above_idx))); }
auto below(NeighborhoodIteratorType const& n) const { return otb::Below<InputV3>(XYZ(n.GetPixel(below_idx))); }
auto left (NeighborhoodIteratorType const& n) const { return otb::Left <InputV3>(XYZ(n.GetPixel(left_idx))); }
auto right(NeighborhoodIteratorType const& n) const { return otb::Right<InputV3>(XYZ(n.GetPixel(right_idx))); }
static constexpr auto center_idx = 4;
auto center(NeighborhoodIteratorType const& n) const { return XYZ(n.GetPixel(center_idx)); }
private:
XYZ_t<InputImageType> XYZ;
int above_idx;
int below_idx;
int left_idx;
int right_idx;
};
/** Specialization that returns pixel means on all four sides. */
template <typename InputImageType>
struct side_extractors_t<InputImageType, true>
{
side_extractors_t(
FloatType x_spacing,
FloatType y_spacing,
FloatType nodata,
InputImageType & image)
: XYZ(image)
, nodata_in{nodata, nodata, nodata}
, all_above_idx{y_spacing > 0 ? std::array<int, 3>{6, 7, 8} : std::array<int, 3>{0, 1, 2}}
, all_below_idx{y_spacing > 0 ? std::array<int, 3>{0, 1, 2} : std::array<int, 3>{6, 7, 8}}
, all_left_idx {x_spacing > 0 ? std::array<int, 3>{0, 3, 6} : std::array<int, 3>{2, 5, 8}}
, all_right_idx{x_spacing > 0 ? std::array<int, 3>{2, 5, 8} : std::array<int, 3>{0, 3, 6}}
{}
using NeighborhoodIteratorType = typename input_types<InputImageType>::NeighborhoodIteratorType;
using InputV3 = typename input_types<InputImageType>::InputV3;
template <typename PixelType>
bool is_invalid0(PixelType const& v) const
{
return
// SARCartesianMeanEstimation will return {0, 0, 0}
(v[0] == 0 && v[1] == 0 && v[2] == 0)
// ossim may produce nan
|| std::isnan(v[0]) || std::isnan(v[1]) || std::isnan(v[2])
#if 0
// other tools may produce nodata;
// but we only support with SARDEMProjection => #if 0
|| v[0] == nodata || v[1] == nodata || v[2] == nodata
#endif
;
}
vnl_vector_fixed<FloatType, 3> mean_side(
NeighborhoodIteratorType const& n,
std::array<int,3> const& all_idx) const
{
unsigned int nb_valids = 0;
vnl_vector_fixed<FloatType, 3> res {0, 0, 0};
for (auto idx : all_idx)
{
auto v = n.GetPixel(idx);
if (! is_invalid0(v))
{
++nb_valids;
res[0] += v[0];
res[1] += v[1];
res[2] += v[2];
}
}
return (nb_valids > 0) ? res / FloatType(nb_valids) : nodata_in;
}
auto above(NeighborhoodIteratorType const& n) const { return otb::Above<InputV3>(mean_side(n, all_above_idx)); }
auto below(NeighborhoodIteratorType const& n) const { return otb::Below<InputV3>(mean_side(n, all_below_idx)); }
auto left (NeighborhoodIteratorType const& n) const { return otb::Left <InputV3>(mean_side(n, all_left_idx )); }
auto right(NeighborhoodIteratorType const& n) const { return otb::Right<InputV3>(mean_side(n, all_right_idx)); }
static constexpr auto center_idx = 4;
auto center(NeighborhoodIteratorType const& n) const { return XYZ(n.GetPixel(center_idx)); }
private:
XYZ_t<InputImageType> XYZ;
vnl_vector_fixed<FloatType, 3> nodata_in;
std::array<int,3> all_above_idx;
std::array<int,3> all_below_idx;
std::array<int,3> all_left_idx;
std::array<int,3> all_right_idx;
};
/** Helper function emulating CTAD in C++14. */
template < bool shall_compute_mean_on_sides, typename InputImageType >
auto make_side_extractors(
FloatType x_spacing,
FloatType y_spacing,
FloatType nodata,
InputImageType & image)
{
// x_spacing > 0 and y_spacing > 0 ==>
// 6 7 8
// 3 4 5
@@ -85,23 +210,40 @@ auto AppSimpleNormalMethod(
// 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));};
return side_extractors_t<InputImageType, shall_compute_mean_on_sides>(
x_spacing, y_spacing, nodata, image);
}
/**
* Returns the lambda that'll be used to produce the normal vector at the
* center of the neighbourhood.
*
* \note As a simplification, this function 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 <bool shall_compute_mean_on_sides,
typename OutputPixelType, typename InputImageType>
inline
auto AppSimpleNormalMethod(
FloatType x_spacing,
FloatType y_spacing,
FloatType nodata,
InputImageType & image)
{
using InputV3 = typename input_types<InputImageType>::InputV3;
using NeighborhoodIteratorType = typename input_types<InputImageType>::NeighborhoodIteratorType;
auto const side_extractors
= make_side_extractors<shall_compute_mean_on_sides>(x_spacing, y_spacing, nodata, image);
auto l0 = otb::SimpleNormalMethodForCartesianPoints<vnl_vector_fixed<FloatType, 3>>(nodata);
auto const nodata_res = otb::repack<OutputPixelType>(nodata, nodata, nodata);
// Need to copy all local variables captured by the lambda
auto lambda = [=](NeighborhoodIteratorType const& in)
{
// Compute n-> at current point
@@ -123,22 +265,18 @@ auto AppSimpleNormalMethod(
;
};
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))
auto const a = side_extractors.above(in);
auto const l = side_extractors.left(in);
auto const r = side_extractors.right(in);
auto const b = side_extractors.below(in);
auto const c = side_extractors.center(in);
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))
{
return nodata_res;
}
@@ -226,6 +364,9 @@ private:
AddParameter(ParameterType_OutputImage, "out", "Image of normal vectors");
SetParameterDescription("out", "Image of normal vectors");
AddParameter(ParameterType_Bool, "mean", "Compute normals from means of cartesian positions on each side");
SetParameterDescription("mean", "Compute normals from means of cartesian positions on each side. By default, use only the positions of the 4 immediate side pixels (right, left, above, below)");
DoInit_NoData();
AddRAMParameter();
@@ -239,6 +380,25 @@ private:
{
}
template <bool does_compute_mean_on_sides>
void RegisterFilter(
FloatType x_spacing, FloatType y_spacing, FloatType nodata,
InputImageType & inputImage)
{
// This function has been made template inorder to avoid branching on pixel
// fetching. Instead the branching is done once when registering the exact
// filter.
auto const lambda = AppSimpleNormalMethod<does_compute_mean_on_sides, 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();
}
void DoExecute() override
{
auto* inputImage = GetParameterImage<InputImageType>("xyz");
@@ -246,21 +406,23 @@ private:
auto const nodata = GetParameterFloat("nodata");
// std::cout << "nodata => " << nodata << std::endl;
bool const does_compute_mean_on_sides = GetParameterInt("mean");
otbLogMacro(Info, <<"Computing normals from "
<< (does_compute_mean_on_sides
? "all pixels on each side" : "a single pixel from each side"));
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);
if (does_compute_mean_on_sides)
RegisterFilter<true>(x_spacing, y_spacing, nodata, *inputImage);
else
RegisterFilter<false>(x_spacing, y_spacing, nodata, *inputImage);
// Call register pipeline to allow streaming and garbage collection
RegisterPipeline();
// RegisterPipeline();
}
};