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

Closes #3 (closed)

app/otbExtractNormalVector.cxx

@@ -89,13 +89,64 @@ auto AppSimpleNormalMethod(
   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;
+
+  // For sides, we don't really care about orders
+  auto const all_above_idx = y_spacing > 0 ? std::array<int, 3>{6, 7, 8} : std::array<int, 3>{0, 1, 2};
+  auto const all_below_idx = y_spacing > 0 ? std::array<int, 3>{0, 1, 2} : std::array<int, 3>{6, 7, 8};
+  auto const all_left_idx  = x_spacing > 0 ? std::array<int, 3>{0, 3, 6} : std::array<int, 3>{2, 5, 8};
+  auto const all_right_idx = x_spacing > 0 ? std::array<int, 3>{2, 5, 8} : std::array<int, 3>{0, 3, 6};
+
   auto const center_idx = 4;
   // TODO: Can we be sure this is enough to guarantee dot(res, center) > 0?
 
+#if 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)));};
+#else
+  auto const is_invalid0 = [nodata](auto v)
+  {
+    // InputV3 const& w = v; // Remove the StrongType<> layer
+    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
+      ;
+  };
+
+  // auto const nodata_in = otb::repack<InputV3>(nodata, nodata, nodata);
+  auto const nodata_in = otb::repack<vnl_vector_fixed<FloatType, 3>>(nodata, nodata, nodata);
+  auto const mean_side  = [is_invalid0, nodata_in, XYZ](NeighborhoodIteratorType const& n, std::array<int, 3> const& all_idx)
+    -> vnl_vector_fixed<FloatType, 3>
+  {
+    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 const above  = [all_above_idx, mean_side, XYZ](NeighborhoodIteratorType const& n) { return otb::Above<InputV3>(mean_side(n, all_above_idx)); };
+  auto const below  = [all_below_idx, mean_side, XYZ](NeighborhoodIteratorType const& n) { return otb::Below<InputV3>(mean_side(n, all_below_idx)); };
+  auto const right  = [all_right_idx, mean_side, XYZ](NeighborhoodIteratorType const& n) { return otb::Right<InputV3>(mean_side(n, all_right_idx)); };
+  auto const left   = [all_left_idx,  mean_side, XYZ](NeighborhoodIteratorType const& n) { return otb::Left <InputV3>(mean_side(n, all_left_idx)); };
+#endif
   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);