Handle no data in SARCalibration

Issue description

Current version of SARCalibration (v9.0) doesn't permit to distinguish nodata pixels from pixels whose physical value after noise removal drops to 0.

"Starting with IPF 2.90, the no-value pixels at image borders are set to 0 value by the IPF. The masking of the borders thus does not require specific processing."

-- https://sentinels.copernicus.eu/documents/247904/2142675/Sentinel-1-masking-no-value-pixels-grd-products-note.pdf/32f11e6f-68b1-4f0a-869b-8d09f80e6788?t=1518545526000

In S1Tiling, a workaround has been implemented to make sure nodata on the sides stays at 0, while pixel noise-corrected to 0 are fixed to a LOWER_SIGNAL_VALUE (= 1e-7 by default). See s1-tiling/s1tiling#87 (closed) and s1-tiling/s1tiling#159 (closed)

Supporting those kinds of nodata situations directly in SARCalibration would simplify S1TIling workflow, but also the workflow of many other users.

Possible solutions

About nodata

Given S1 products, we could just continue to assume nodata to always be 0. But... may be, it'd be better to have a -nodata parameter defaulted to 0. I don't know.

About LOWER_SIGNAL_VALUE

IMO, we should just force the final output to LOWER_SIGNAL_VALUE option, when set (it has to be optional to keep the previous behaviour as the default behaviour), and when -removenoise is True.

Possible implementation

SarRadiometricCalibrationFunction<TInputImage, TCoordRep>::EvaluateAtIndex(const IndexType& index) const
{
...
  const std::complex<float> pVal           = this->GetInputImage()->GetPixel(index);
  // as long as the nodata value is 0, it will be simple with complex numbers
  // that's why I'm not sure we need to handle cases like nodata = -32768 ...
  if (pVal == m_nodata) return m_nodata;
  const RealType            digitalNumber2 = (pVal.real() * pVal.real()) + (pVal.imag() * pVal.imag());
  // why taking the sqrt of the previous value to multiply it by itself? ...
  // see a future Issue I'll open.

...
  if (m_ApplyLookupDataCorrection)
  {
    if (m_EnableNoise && m_NoiseLut)
    {
      sigma = otb::quiet::max(0., sigma - m_NoiseLut->GetValue(index[0], index[1]));
      // See next section about NaN's regarding otb::quiet::max.

      // Given lower_signal_value option is set, and result would have been 0.0
      if (sigma == 0. && m_using_lower_signal_value)
        return m_lower_signal_value;
      // no need to handle NaN values as further computations will continue to
      // return NaNs... But we could also guard the function for this case.
    }

    RealType lutVal = static_cast<RealType>(m_Lut->GetValue(index[0], index[1]));

    sigma /= lutVal * lutVal;
  }
...

About NaN values

Given peculiar handling of NaN numbers in C and C++ (1), we tried to use BandMath to replace all 0's from an original S1 image with NaNs, thanks to exp='im1b1==im1b1?im1b1:0./0. (See #1798 (comment 110682))

SARCalibration without noise removal keeps the NaN on the sides (which is what we were looking for). However as soon as noiseremoval is requested, NaN values are replaced with 0's.

This is because of the following:

// otbSarRadiometricCalibrationFunction.hxx
// SarRadiometricCalibrationFunction<>::EvaluateAtIndex {...
    if (m_EnableNoise && m_NoiseLut)
    {
      sigma = std::max(0., sigma - m_NoiseLut->GetValue(index[0], index[1]));
    }

Indeed with current implementation of std::max in GNU libstdc++, std::max(whatever, nan) will always return whatever. On the contrary, std::max(nan, whatever) will always return nan.

    // current implementation of std::max in libstdc++
    max(const _Tp& __a, const _Tp& __b)
    {
      if (__a < __b) // always false when a or b is NaN
        return __b;
      return __a;
    }

Also, it seems the behaviour of std::max(nan, ...) is undefined ("UB"). In consequence, we have to take control over max() definition to guarantee a defined behaviour when NaN values are passed.

=> recommendation

See: https://godbolt.org/z/633Y5ojT7

1. Define a flavour of max without UB on NaN values

// The following is strongly inspired from
// Walter E. Brown - Correctly calculating min, max and more - Meeting C++ online
// https://www.youtube.com/watch?v=V1U6AXEgNbE


// PS: I've used lambdas instead of functions to simplify the definitions.
// Inline (!!!) functions would be perfect.

namespace quiet {
auto out_of_order = [](double a, double b) {
  // std::isless has a guaranteed non signalling behaviour when NaN values are passed
  return std::isless(b, a);
};
auto max = [](double a, double b) {
  return out_of_order(a, b) ? a : b;
};
}

2. Update SARCalibration code

// Replace
      sigma = std::max(0., sigma - m_NoiseLut->GetValue(index[0], index[1]));
// With
      sigma = otb::quiet::max(0., sigma - m_NoiseLut->GetValue(index[0], index[1]));

Note: we shall not change parameters orders from the moment max() is "correctly" implemented.

Or...

Or... we could simply guard the function. Still, this NaN related issue may appear in other OTB filters and applications -- we have never tried to correctly handle NaN values in OTB as an alternate nodata value.

So. Instead of

if (pVal == m_nodata) return m_nodata;

we could guard the function with:

if (std::isnan(pVal) || pVal == m_nodata) return pVal;

(1) What's so special about NaN

• Arithmetic functions on nan's will return nan
• comparisons on NaN will always return false, but they may also raise FE_INVALID; this is why std::isless can be preferred.