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 
31 namespace 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 
38 template <typename S, typename T>
39 class Interpolator : public std::map<S, T> {
40 public:
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)
mCache(cache)54         : mCache(cache)
55         , mFirstSlope(0)
56         , mLastSlope(0) {
57         setInterpolatorType(interpolatorType);
58     }
59 
first()60     std::pair<S, T> first() const {
61         return *this->begin();
62     }
63 
last()64     std::pair<S, T> last() const {
65         return *this->rbegin();
66     }
67 
68     // find the corresponding Y point from a X point.
findY(S x)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 
getInterpolatorType()199     InterpolatorType getInterpolatorType() const {
200         return mInterpolatorType;
201     }
202 
setInterpolatorType(InterpolatorType interpolatorType)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 
getFirstSlope()218     T getFirstSlope() const {
219         return mFirstSlope;
220     }
221 
setFirstSlope(T slope)222     void setFirstSlope(T slope) {
223         mFirstSlope = slope;
224     }
225 
getLastSlope()226     T getLastSlope() const {
227         return mLastSlope;
228     }
229 
setLastSlope(T slope)230     void setLastSlope(T slope) {
231         mLastSlope = slope;
232     }
233 
clearCache()234     void clearCache() {
235         mMemo.clear();
236     }
237 
writeToParcel(Parcel * parcel)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 
readFromParcel(const Parcel & parcel)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 
toString()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 
315 private:
constrainSlope(S slope,S maxSlope)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