1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
|
#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;
}
|