Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • s1-tiling/normlim_sigma0
1 result
Show changes
Commits on Source (4)
Hide whitespace changes
Inline Side-by-side
app/otbExtractNormalVector.cxx
View file @ 6196399e
......@@ -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();
}
};
......