innards.impl.h

Go to the documentation of this file.

#ifndef SOPT_WAVELETS_INNARDS_H
#define SOPT_WAVELETS_INNARDS_H
 
#include "sopt/config.h"
#include <algorithm> // for std::copy<>
#include <Eigen/Core>
 
// Function inside anonymouns namespace won't appear in library
namespace sopt::wavelets {
namespace {
template <typename T0>
auto copy(Eigen::ArrayBase<T0> const &a) ->
    typename std::remove_const<typename std::remove_reference<decltype(a.eval())>::type>::type {
  return a.eval();
}
 
template <typename T0, typename T2>
typename T0::Scalar periodic_scalar_product(Eigen::ArrayBase<T0> const &a,
                                            Eigen::ArrayBase<T2> const &b,
                                            typename T0::Index offset) {
  auto const Na = static_cast<typename T0::Index>(a.size());
  auto const Nb = static_cast<typename T0::Index>(b.size());
  // offset in [0, a.size()[
  offset %= Na;
  if (offset < 0) offset += Na;
  // Simple case, just do it
  if (Na > Nb + offset) {
    return (a.segment(offset, Nb) * b).sum();
  }
  // Wrap around, do it, but carefully
  typename T0::Scalar result(0);
  for (typename T0::Index i(0), j(offset); i < Nb; ++i, ++j) result += a(j % Na) * b(i);
  return result;
}
 
template <typename T0, typename T1, typename T2>
void convolve(Eigen::ArrayBase<T0> &result, Eigen::ArrayBase<T1> const &signal,
              Eigen::ArrayBase<T2> const &filter) {
  assert(result.size() == signal.size());
  for (typename T0::Index i(0); i < t_int(signal.size()); ++i)
    result(i) = periodic_scalar_product(signal, filter, i);
}
template <typename T0, typename T1, typename T2>
void convolve(Eigen::VectorBlock<T0> &&result, Eigen::ArrayBase<T1> const &signal,
              Eigen::ArrayBase<T2> const &filter) {
  return convolve(result, signal, filter);
}
 
template <typename T0, typename T1, typename T2>
void down_convolve(Eigen::ArrayBase<T0> &result, Eigen::ArrayBase<T1> const &signal,
                   Eigen::ArrayBase<T2> const &filter) {
  assert(result.size() * 2 <= signal.size());
  if (signal.rows() == 1) {
#ifdef SOPT_OPENMP
#pragma omp parallel for
#endif
    for (typename T0::Index i = 0; i < result.size(); ++i)
      result(i) = periodic_scalar_product(signal.transpose(), filter, 2 * i);
  } else {
#ifdef SOPT_OPENMP
#pragma omp parallel for
#endif
    for (typename T0::Index i = 0; i < result.size(); ++i)
      result(i) = periodic_scalar_product(signal, filter, 2 * i);
  }
}
template <typename T0, typename T1, typename T2>
void down_convolve(Eigen::VectorBlock<T0> &&result, Eigen::ArrayBase<T1> const &signal,
                   Eigen::ArrayBase<T2> const &filter) {
  down_convolve(result, signal, filter);
}
 
template <typename T0, typename T1, typename T2, typename T3, typename T4>
void convolve_sum(Eigen::ArrayBase<T0> &result, Eigen::ArrayBase<T1> const &low_pass_signal,
                  Eigen::ArrayBase<T2> const &low_pass,
                  Eigen::ArrayBase<T3> const &high_pass_signal,
                  Eigen::ArrayBase<T4> const &high_pass) {
  assert(result.size() == low_pass_signal.size());
  assert(result.size() == high_pass_signal.size());
  assert(low_pass.size() == high_pass.size());
  static_assert(std::is_signed<typename T0::Index>::value, "loffset, hoffset expect signed values");
  auto const loffset = 1 - static_cast<typename T0::Index>(low_pass.size());
  auto const hoffset = 1 - static_cast<typename T0::Index>(high_pass.size());
  for (typename T0::Index i(0); i < result.size(); ++i) {
    result(i) = periodic_scalar_product(low_pass_signal, low_pass, i + loffset) +
                periodic_scalar_product(high_pass_signal, high_pass, i + hoffset);
  }
}
 
template <typename T0, typename T1, typename T2, typename T3, typename T4, typename T5>
void up_convolve_sum(Eigen::ArrayBase<T0> &result, Eigen::ArrayBase<T1> const &coeffs,
                     Eigen::ArrayBase<T2> const &low_even, Eigen::ArrayBase<T3> const &low_odd,
                     Eigen::ArrayBase<T4> const &high_even, Eigen::ArrayBase<T5> const &high_odd) {
  assert(result.size() <= coeffs.size());
  assert(result.size() % 2 == 0);
  assert(low_even.size() == high_even.size());
  assert(low_odd.size() == high_odd.size());
  auto const Nlow = (coeffs.size() + 1) / 2;
  auto const Nhigh = coeffs.size() / 2;
  auto const size = low_even.size() + low_odd.size();
  auto const is_even = size % 2 == 0;
  auto const even_offset = (1 - size) / 2;
  auto const odd_offset = (1 - size) / 2 + (is_even ? 0 : 1);
  auto const index_place_even = (is_even ? 1 : 0);
  auto const index_place_odd = (is_even ? 0 : 1);
#ifdef SOPT_OPENMP
#pragma omp parallel for
#endif
  for (typename T0::Index i = 0; i < result.size() / 2; i++) {
    result(2 * i + index_place_even) =
        periodic_scalar_product(coeffs.head(Nlow), low_even, i + even_offset) +
        periodic_scalar_product(coeffs.tail(Nhigh), high_even, i + even_offset);
    result(2 * i + index_place_odd) =
        periodic_scalar_product(coeffs.head(Nlow), low_odd, i + odd_offset) +
        periodic_scalar_product(coeffs.tail(Nhigh), high_odd, i + odd_offset);
  }
}
}  // namespace
} // namespace sopt::wavelets
#endif