1/*
2 * Copyright 2017 The Android Open Source Project
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 *      http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17#ifndef ANDROID_INTERPOLATOR_H
18#define ANDROID_INTERPOLATOR_H
19
20#include <map>
21#include <sstream>
22#include <unordered_map>
23
24#include <binder/Parcel.h>
25#include <utils/RefBase.h>
26
27#pragma push_macro("LOG_TAG")
28#undef LOG_TAG
29#define LOG_TAG "Interpolator"
30
31namespace android {
32
33/*
34 * A general purpose spline interpolator class which takes a set of points
35 * and performs interpolation.  This is used for the VolumeShaper class.
36 */
37
38template <typename S, typename T>
39class Interpolator : public std::map<S, T> {
40public:
41    // Polynomial spline interpolators
42    // Extend only at the end of enum, as this must match order in VolumeShapers.java.
43    enum InterpolatorType : int32_t {
44        INTERPOLATOR_TYPE_STEP,   // Not continuous
45        INTERPOLATOR_TYPE_LINEAR, // C0
46        INTERPOLATOR_TYPE_CUBIC,  // C1
47        INTERPOLATOR_TYPE_CUBIC_MONOTONIC, // C1 (to provide locally monotonic curves)
48        // INTERPOLATOR_TYPE_CUBIC_C2, // TODO - requires global computation / cache
49    };
50
51    explicit Interpolator(
52            InterpolatorType interpolatorType = INTERPOLATOR_TYPE_LINEAR,
53            bool cache = true)
54        : mCache(cache)
55        , mFirstSlope(0)
56        , mLastSlope(0) {
57        setInterpolatorType(interpolatorType);
58    }
59
60    std::pair<S, T> first() const {
61        return *this->begin();
62    }
63
64    std::pair<S, T> last() const {
65        return *this->rbegin();
66    }
67
68    // find the corresponding Y point from a X point.
69    T findY(S x) { // logically const, but modifies cache
70        auto high = this->lower_bound(x);
71        // greater than last point
72        if (high == this->end()) {
73            return this->rbegin()->second;
74        }
75        // at or before first point
76        if (high == this->begin()) {
77            return high->second;
78        }
79        // go lower.
80        auto low = high;
81        --low;
82
83        // now that we have two adjacent points:
84        switch (mInterpolatorType) {
85        case INTERPOLATOR_TYPE_STEP:
86            return high->first == x ? high->second : low->second;
87        case INTERPOLATOR_TYPE_LINEAR:
88            return ((high->first - x) * low->second + (x - low->first) * high->second)
89                    / (high->first - low->first);
90        case INTERPOLATOR_TYPE_CUBIC:
91        case INTERPOLATOR_TYPE_CUBIC_MONOTONIC:
92        default: {
93            // See https://en.wikipedia.org/wiki/Cubic_Hermite_spline
94
95            const S interval =  high->first - low->first;
96
97            // check to see if we've cached the polynomial coefficients
98            if (mMemo.count(low->first) != 0) {
99                const S t = (x - low->first) / interval;
100                const S t2 = t * t;
101                const auto &memo = mMemo[low->first];
102                return low->second + std::get<0>(memo) * t
103                        + (std::get<1>(memo) + std::get<2>(memo) * t) * t2;
104            }
105
106            // find the neighboring points (low2 < low < high < high2)
107            auto low2 = this->end();
108            if (low != this->begin()) {
109                low2 = low;
110                --low2; // decrementing this->begin() is undefined
111            }
112            auto high2 = high;
113            ++high2;
114
115            // you could have catmullRom with monotonic or
116            // non catmullRom (finite difference) with regular cubic;
117            // the choices here minimize computation.
118            bool monotonic, catmullRom;
119            if (mInterpolatorType == INTERPOLATOR_TYPE_CUBIC_MONOTONIC) {
120                monotonic = true;
121                catmullRom = false;
122            } else {
123                monotonic = false;
124                catmullRom = true;
125            }
126
127            // secants are only needed for finite difference splines or
128            // monotonic computation.
129            // we use lazy computation here - if we precompute in
130            // a single pass, duplicate secant computations may be avoided.
131            S sec, sec0, sec1;
132            if (!catmullRom || monotonic) {
133                sec = (high->second - low->second) / interval;
134                sec0 = low2 != this->end()
135                        ? (low->second - low2->second) / (low->first - low2->first)
136                        : mFirstSlope;
137                sec1 = high2 != this->end()
138                        ? (high2->second - high->second) / (high2->first - high->first)
139                        : mLastSlope;
140            }
141
142            // compute the tangent slopes at the control points
143            S m0, m1;
144            if (catmullRom) {
145                // Catmull-Rom spline
146                m0 = low2 != this->end()
147                        ? (high->second - low2->second) / (high->first - low2->first)
148                        : mFirstSlope;
149
150                m1 = high2 != this->end()
151                        ? (high2->second - low->second) / (high2->first - low->first)
152                        : mLastSlope;
153            } else {
154                // finite difference spline
155                m0 = (sec0 + sec) * 0.5f;
156                m1 = (sec1 + sec) * 0.5f;
157            }
158
159            if (monotonic) {
160                // https://en.wikipedia.org/wiki/Monotone_cubic_interpolation
161                // A sufficient condition for Fritsch–Carlson monotonicity is constraining
162                // (1) the normalized slopes to be within the circle of radius 3, or
163                // (2) the normalized slopes to be within the square of radius 3.
164                // Condition (2) is more generous and easier to compute.
165                const S maxSlope = 3 * sec;
166                m0 = constrainSlope(m0, maxSlope);
167                m1 = constrainSlope(m1, maxSlope);
168
169                m0 = constrainSlope(m0, 3 * sec0);
170                m1 = constrainSlope(m1, 3 * sec1);
171            }
172
173            const S t = (x - low->first) / interval;
174            const S t2 = t * t;
175            if (mCache) {
176                // convert to cubic polynomial coefficients and compute
177                m0 *= interval;
178                m1 *= interval;
179                const T dy = high->second - low->second;
180                const S c0 = low->second;
181                const S c1 = m0;
182                const S c2 = 3 * dy - 2 * m0 - m1;
183                const S c3 = m0 + m1 - 2 * dy;
184                mMemo[low->first] = std::make_tuple(c1, c2, c3);
185                return c0 + c1 * t + (c2 + c3 * t) * t2;
186            } else {
187                // classic Hermite interpolation
188                const S t3 = t2 * t;
189                const S h00 =  2 * t3 - 3 * t2     + 1;
190                const S h10 =      t3 - 2 * t2 + t    ;
191                const S h01 = -2 * t3 + 3 * t2        ;
192                const S h11 =      t3     - t2        ;
193                return h00 * low->second + (h10 * m0 + h11 * m1) * interval + h01 * high->second;
194            }
195        } // default
196        }
197    }
198
199    InterpolatorType getInterpolatorType() const {
200        return mInterpolatorType;
201    }
202
203    status_t setInterpolatorType(InterpolatorType interpolatorType) {
204        switch (interpolatorType) {
205        case INTERPOLATOR_TYPE_STEP:   // Not continuous
206        case INTERPOLATOR_TYPE_LINEAR: // C0
207        case INTERPOLATOR_TYPE_CUBIC:  // C1
208        case INTERPOLATOR_TYPE_CUBIC_MONOTONIC: // C1 + other constraints
209        // case INTERPOLATOR_TYPE_CUBIC_C2:
210            mInterpolatorType = interpolatorType;
211            return NO_ERROR;
212        default:
213            ALOGE("invalid interpolatorType: %d", interpolatorType);
214            return BAD_VALUE;
215        }
216    }
217
218    T getFirstSlope() const {
219        return mFirstSlope;
220    }
221
222    void setFirstSlope(T slope) {
223        mFirstSlope = slope;
224    }
225
226    T getLastSlope() const {
227        return mLastSlope;
228    }
229
230    void setLastSlope(T slope) {
231        mLastSlope = slope;
232    }
233
234    void clearCache() {
235        mMemo.clear();
236    }
237
238    status_t writeToParcel(Parcel *parcel) const {
239        if (parcel == nullptr) {
240            return BAD_VALUE;
241        }
242        status_t res = parcel->writeInt32(mInterpolatorType)
243                ?: parcel->writeFloat(mFirstSlope)
244                ?: parcel->writeFloat(mLastSlope)
245                ?: parcel->writeUint32((uint32_t)this->size()); // silent truncation
246        if (res != NO_ERROR) {
247            return res;
248        }
249        for (const auto &pt : *this) {
250            res = parcel->writeFloat(pt.first)
251                    ?: parcel->writeFloat(pt.second);
252            if (res != NO_ERROR) {
253                return res;
254            }
255        }
256        return NO_ERROR;
257    }
258
259    status_t readFromParcel(const Parcel &parcel) {
260        this->clear();
261        int32_t type;
262        uint32_t size;
263        status_t res = parcel.readInt32(&type)
264                        ?: parcel.readFloat(&mFirstSlope)
265                        ?: parcel.readFloat(&mLastSlope)
266                        ?: parcel.readUint32(&size)
267                        ?: setInterpolatorType((InterpolatorType)type);
268        if (res != NO_ERROR) {
269            return res;
270        }
271        // Note: We don't need to check size is within some bounds as
272        // the Parcel read will fail if size is incorrectly specified too large.
273        float lastx;
274        for (uint32_t i = 0; i < size; ++i) {
275            float x, y;
276            res = parcel.readFloat(&x)
277                    ?: parcel.readFloat(&y);
278            if (res != NO_ERROR) {
279                return res;
280            }
281            if (i > 0 && !(x > lastx) /* handle nan */
282                    || y != y /* handle nan */) {
283                // This is a std::map object which imposes sorted order
284                // automatically on emplace.
285                // Nevertheless for reading from a Parcel,
286                // we require that the points be specified monotonic in x.
287                return BAD_VALUE;
288            }
289            this->emplace(x, y);
290            lastx = x;
291        }
292        return NO_ERROR;
293    }
294
295    std::string toString() const {
296        std::stringstream ss;
297        ss << "Interpolator{mInterpolatorType=" << static_cast<int32_t>(mInterpolatorType);
298        ss << ", mFirstSlope=" << mFirstSlope;
299        ss << ", mLastSlope=" << mLastSlope;
300        ss << ", {";
301        bool first = true;
302        for (const auto &pt : *this) {
303            if (first) {
304                first = false;
305                ss << "{";
306            } else {
307                ss << ", {";
308            }
309            ss << pt.first << ", " << pt.second << "}";
310        }
311        ss << "}}";
312        return ss.str();
313    }
314
315private:
316    static S constrainSlope(S slope, S maxSlope) {
317        if (maxSlope > 0) {
318            slope = std::min(slope, maxSlope);
319            slope = std::max(slope, S(0)); // not globally monotonic
320        } else {
321            slope = std::max(slope, maxSlope);
322            slope = std::min(slope, S(0)); // not globally monotonic
323        }
324        return slope;
325    }
326
327    InterpolatorType mInterpolatorType;
328    bool mCache; // whether we cache spline coefficient computation
329
330    // for cubic interpolation, the boundary conditions in slope.
331    S mFirstSlope;
332    S mLastSlope;
333
334    // spline cubic polynomial coefficient cache
335    std::unordered_map<S, std::tuple<S /* c1 */, S /* c2 */, S /* c3 */>> mMemo;
336}; // Interpolator
337
338} // namespace android
339
340#pragma pop_macro("LOG_TAG")
341
342#endif // ANDROID_INTERPOLATOR_H
343