KernelSmoother.h

Go to the documentation of this file.

// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
// SPDX-License-Identifier: GPL-3.0-or-later
 
#ifndef DP3_DDECAL_KERNEL_SMOOTHER_H_
#define DP3_DDECAL_KERNEL_SMOOTHER_H_
 
#include <cmath>
#include <complex>
#include <stdexcept>
#include <vector>
#include <limits>
 
namespace dp3 {
namespace ddecal {
 
template <typename DataType, typename NumType>
class KernelSmoother {
 public:
  enum KernelType {
    RectangularKernel,
    TriangularKernel,
    GaussianKernel,
    EpanechnikovKernel
  };
 
  KernelSmoother(const std::vector<NumType>& frequencies, KernelType kernelType,
                 NumType kernelBandwidth, NumType bandwidthRefFrequency,
                 NumType spectralExponent, bool kernel_truncation)
      : _frequencies(frequencies),
        _scratch(frequencies.size()),
        _kernelType(kernelType),
        _bandwidth(kernelBandwidth),
        _bandwidthRefFrequency(bandwidthRefFrequency),
        _spectralExponent(spectralExponent),
        _truncate(kernel_truncation) {}
 
  NumType Kernel(NumType distance) const {
    NumType x = distance / _bandwidth;
    switch (_kernelType) {
      case RectangularKernel:
      default:
        if (x < NumType(-1.0) || x > NumType(1.0))
          return NumType(0.0);
        else
          return NumType(0.5);
      case TriangularKernel:
        if (x < NumType(-1.0) || x > NumType(1.0))
          return NumType(0.0);
        else
          return x >= NumType(0.0) ? (NumType(1.0) - x) : (NumType(1.0) + x);
      case GaussianKernel:
        if (_truncate && (x < NumType(-1.0) || x > NumType(1.0)))
          return 0.0;
        else
          return std::exp(-x * x * NumType(9.0));
      case EpanechnikovKernel:
        if (_truncate && (x < NumType(-1.0) || x > NumType(1.0))) {
          return 0.0;
        } else {
          x = NumType(1.0) - x;
          return (NumType(3.0) / NumType(4.0)) * x * x;
        }
    }
  }
 
  void Smooth(DataType* data, const NumType* weight, NumType kernelSizeFactor) {
    size_t n = _frequencies.size();
 
    size_t bandLeft = 0;
    const size_t lbound =
        std::lower_bound(_frequencies.begin(), _frequencies.end(),
                         _frequencies[0] + _bandwidth * 0.5) -
        _frequencies.begin();
    size_t bandRight = std::min(lbound + 1u, n);
    for (size_t i = 0; i != n; ++i) {
      size_t start = 0;
      size_t end = n;
      const NumType localBandwidth =
          _bandwidthRefFrequency == 0.0
              ? _bandwidth
              : _bandwidth * _frequencies[i] / _bandwidthRefFrequency;
      if (_truncate) {
        while (_frequencies[bandLeft] < _frequencies[i] - localBandwidth * 0.5)
          ++bandLeft;
        while (bandRight != n &&
               _frequencies[bandRight] < _frequencies[i] + localBandwidth * 0.5)
          ++bandRight;
 
        start = bandLeft > 0 ? bandLeft - 1 : 0;
        end = bandRight < n ? bandRight + 1 : n;
      }
 
      DataType sum(0.0);
      NumType weightSum(0.0);
      const NumType frequencyCorrection =
          std::pow(localBandwidth / _bandwidth, _spectralExponent);
      const NumType kernelCorrection = frequencyCorrection * kernelSizeFactor;
      for (size_t j = start; j != end; ++j) {
        NumType distance = _frequencies[i] - _frequencies[j];
        double w = Kernel(distance * kernelCorrection) * weight[j];
        sum += data[j] * w;
        weightSum += w;
      }
      if (weightSum == 0.0)
        _scratch[i] = quiet_NaN(_scratch[i]);
      else
        _scratch[i] = sum / weightSum;
    }
    std::copy(_scratch.begin(), _scratch.end(), data);
  }
 
 private:
  double quiet_NaN(double) { return std::numeric_limits<double>::quiet_NaN(); }
 
  std::complex<double> quiet_NaN(std::complex<double>) {
    return std::complex<double>(std::numeric_limits<double>::quiet_NaN(),
                                std::numeric_limits<double>::quiet_NaN());
  }
 
  std::vector<NumType> _frequencies;
  std::vector<DataType> _scratch;
  enum KernelType _kernelType;
  NumType _bandwidth;
  NumType _bandwidthRefFrequency;
  NumType _spectralExponent;
  bool _truncate = true;
};
 
}  // namespace ddecal
}  // namespace dp3
 
#endif