osb/source/core/StarMultiTable.hpp

#pragma once

#include "StarMultiArrayInterpolator.hpp"

namespace Star {

// Provides a method for storing, retrieving, and interpolating uneven
// n-variate data.  Access times involve a binary search over the domain of
// each dimension, so is O(log(n)*m) where n is the size of the largest
// dimension, and m is the table_rank.
template <typename ElementT, typename PositionT, size_t RankN>
class MultiTable {
public:
  typedef ElementT Element;
  typedef PositionT Position;
  static size_t const Rank = RankN;

  typedef Star::MultiArray<ElementT, RankN> MultiArray;

  typedef Star::MultiArrayInterpolator2<MultiArray, Position> Interpolator2;
  typedef Star::MultiArrayInterpolator4<MultiArray, Position> Interpolator4;
  typedef Star::MultiArrayPiecewiseInterpolator<MultiArray, Position> PiecewiseInterpolator;

  typedef Array<Position, Rank> PositionArray;
  typedef Array<Position, 2> WeightArray2;
  typedef Array<Position, 4> WeightArray4;
  typedef typename MultiArray::SizeArray SizeArray;
  typedef typename MultiArray::IndexArray IndexArray;
  typedef List<Position> Range;
  typedef Array<Range, Rank> RangeArray;

  typedef std::function<WeightArray2(Position)> WeightFunction2;
  typedef std::function<WeightArray4(Position)> WeightFunction4;
  typedef std::function<Element(PositionArray const&)> InterpolateFunction;

  MultiTable() : m_interpolationMode(InterpolationMode::Linear), m_boundMode(BoundMode::Clamp) {}

  // Set input ranges on a particular dimension.  Will resize underlying storage
  // to fit range.
  void setRange(std::size_t dim, Range const& range) {
    SizeArray sizes = m_array.size();
    sizes[dim] = range.size();
    m_array.resize(sizes);

    m_ranges[dim] = range;
  }

  void setRanges(RangeArray const& ranges) {
    SizeArray arraySize;

    for (size_t dim = 0; dim < Rank; ++dim) {
      arraySize[dim] = ranges[dim].size();
      m_ranges[dim] = ranges[dim];
    }

    m_array.resize(arraySize);
  }

  // Set array element based on index.
  void set(IndexArray const& index, Element const& element) {
    m_array.set(index, element);
  }

  // Get array element based on index.
  Element const& get(IndexArray const& index) const {
    return m_array(index);
  }

  MultiArray const& array() const {
    return m_array;
  }

  MultiArray& array() {
    return m_array;
  }

  InterpolationMode interpolationMode() const {
    return m_interpolationMode;
  }

  void setInterpolationMode(InterpolationMode interpolationMode) {
    m_interpolationMode = interpolationMode;
  }

  BoundMode boundMode() const {
    return m_boundMode;
  }

  void setBoundMode(BoundMode boundMode) {
    m_boundMode = boundMode;
  }

  Element interpolate(PositionArray const& coord) const {
    if (m_interpolationMode == InterpolationMode::HalfStep) {
      PiecewiseInterpolator piecewiseInterpolator(StepWeightOperator<Position>(), m_boundMode);
      return piecewiseInterpolator.interpolate(m_array, toIndexSpace(coord));

    } else if (m_interpolationMode == InterpolationMode::Linear) {
      Interpolator2 interpolator2(LinearWeightOperator<Position>(), m_boundMode);
      return interpolator2.interpolate(m_array, toIndexSpace(coord));

    } else if (m_interpolationMode == InterpolationMode::Cubic) {
      // MultiTable uses CubicWeights with linear extrapolation (not
      // configurable atm)
      Interpolator4 interpolator4(Cubic4WeightOperator<Position>(true), m_boundMode);
      return interpolator4.interpolate(m_array, toIndexSpace(coord));

    } else {
      throw MathException("Unsupported interpolation type in MultiTable::interpolate");
    }
  }

  // Synonym for inteprolate
  Element operator()(PositionArray const& coord) const {
    return interpolate(coord);
  }

  // op should take a PositionArray parameter and return an element.
  template <typename OpType>
  void eval(OpType op) {
    m_array.forEach(EvalWrapper<OpType>(op, *this));
  }

private:
  template <typename Coordinate>
  inline PositionArray toIndexSpace(Coordinate const& coord) const {
    PositionArray indexCoord;
    for (size_t i = 0; i < Rank; ++i)
      indexCoord[i] = inverseLinearInterpolateLower(m_ranges[i].begin(), m_ranges[i].end(), coord[i]);
    return indexCoord;
  }

  template <typename OpType>
  struct EvalWrapper {
    EvalWrapper(OpType& o, MultiTable const& t)
      : op(o), table(t) {}

    template <typename IndexArray>
    void operator()(IndexArray const& indexArray, Element& element) {
      PositionArray rangeArray;
      for (size_t i = 0; i < Rank; ++i)
        rangeArray[i] = table.m_ranges[i][indexArray[i]];

      element = op(rangeArray);
    }

    OpType& op;
    MultiTable const& table;
  };

  RangeArray m_ranges;
  MultiArray m_array;
  InterpolationMode m_interpolationMode;
  BoundMode m_boundMode;
};

typedef MultiTable<float, float, 2> MultiTable2F;
typedef MultiTable<double, double, 2> MultiTable2D;

typedef MultiTable<float, float, 3> MultiTable3F;
typedef MultiTable<double, double, 3> MultiTable3D;

typedef MultiTable<float, float, 4> MultiTable4F;
typedef MultiTable<double, double, 4> MultiTable4D;

}