include/media/Interpolator.h

/*
 * Copyright 2017 The Android Open Source Project
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef ANDROID_INTERPOLATOR_H
#define ANDROID_INTERPOLATOR_H

#include <map>
#include <sstream>
#include <unordered_map>

#include <binder/Parcel.h>
#include <utils/RefBase.h>

#pragma push_macro("LOG_TAG")
#undef LOG_TAG
#define LOG_TAG "Interpolator"

namespace android {

/*
 * A general purpose spline interpolator class which takes a set of points
 * and performs interpolation.  This is used for the VolumeShaper class.
 */

template <typename S, typename T>
class Interpolator : public std::map<S, T> {
public:
    // Polynomial spline interpolators
    // Extend only at the end of enum, as this must match order in VolumeShapers.java.
    enum InterpolatorType : int32_t {
        INTERPOLATOR_TYPE_STEP,   // Not continuous
        INTERPOLATOR_TYPE_LINEAR, // C0
        INTERPOLATOR_TYPE_CUBIC,  // C1
        INTERPOLATOR_TYPE_CUBIC_MONOTONIC, // C1 (to provide locally monotonic curves)
        // INTERPOLATOR_TYPE_CUBIC_C2, // TODO - requires global computation / cache
    };

    explicit Interpolator(
            InterpolatorType interpolatorType = INTERPOLATOR_TYPE_LINEAR,
            bool cache = true)
        : mCache(cache)
        , mFirstSlope(0)
        , mLastSlope(0) {
        setInterpolatorType(interpolatorType);
    }

    std::pair<S, T> first() const {
        return *this->begin();
    }

    std::pair<S, T> last() const {
        return *this->rbegin();
    }

    // find the corresponding Y point from a X point.
    T findY(S x) { // logically const, but modifies cache
        auto high = this->lower_bound(x);
        // greater than last point
        if (high == this->end()) {
            return this->rbegin()->second;
        }
        // at or before first point
        if (high == this->begin()) {
            return high->second;
        }
        // go lower.
        auto low = high;
        --low;

        // now that we have two adjacent points:
        switch (mInterpolatorType) {
        case INTERPOLATOR_TYPE_STEP:
            return high->first == x ? high->second : low->second;
        case INTERPOLATOR_TYPE_LINEAR:
            return ((high->first - x) * low->second + (x - low->first) * high->second)
                    / (high->first - low->first);
        case INTERPOLATOR_TYPE_CUBIC:
        case INTERPOLATOR_TYPE_CUBIC_MONOTONIC:
        default: {
            // See https://en.wikipedia.org/wiki/Cubic_Hermite_spline

            const S interval =  high->first - low->first;

            // check to see if we've cached the polynomial coefficients
            if (mMemo.count(low->first) != 0) {
                const S t = (x - low->first) / interval;
                const S t2 = t * t;
                const auto &memo = mMemo[low->first];
                return low->second + std::get<0>(memo) * t
                        + (std::get<1>(memo) + std::get<2>(memo) * t) * t2;
            }

            // find the neighboring points (low2 < low < high < high2)
            auto low2 = this->end();
            if (low != this->begin()) {
                low2 = low;
                --low2; // decrementing this->begin() is undefined
            }
            auto high2 = high;
            ++high2;

            // you could have catmullRom with monotonic or
            // non catmullRom (finite difference) with regular cubic;
            // the choices here minimize computation.
            bool monotonic, catmullRom;
            if (mInterpolatorType == INTERPOLATOR_TYPE_CUBIC_MONOTONIC) {
                monotonic = true;
                catmullRom = false;
            } else {
                monotonic = false;
                catmullRom = true;
            }

            // secants are only needed for finite difference splines or
            // monotonic computation.
            // we use lazy computation here - if we precompute in
            // a single pass, duplicate secant computations may be avoided.
            S sec, sec0, sec1;
            if (!catmullRom || monotonic) {
                sec = (high->second - low->second) / interval;
                sec0 = low2 != this->end()
                        ? (low->second - low2->second) / (low->first - low2->first)
                        : mFirstSlope;
                sec1 = high2 != this->end()
                        ? (high2->second - high->second) / (high2->first - high->first)
                        : mLastSlope;
            }

            // compute the tangent slopes at the control points
            S m0, m1;
            if (catmullRom) {
                // Catmull-Rom spline
                m0 = low2 != this->end()
                        ? (high->second - low2->second) / (high->first - low2->first)
                        : mFirstSlope;

                m1 = high2 != this->end()
                        ? (high2->second - low->second) / (high2->first - low->first)
                        : mLastSlope;
            } else {
                // finite difference spline
                m0 = (sec0 + sec) * 0.5f;
                m1 = (sec1 + sec) * 0.5f;
            }

            if (monotonic) {
                // https://en.wikipedia.org/wiki/Monotone_cubic_interpolation
                // A sufficient condition for Fritsch–Carlson monotonicity is constraining
                // (1) the normalized slopes to be within the circle of radius 3, or
                // (2) the normalized slopes to be within the square of radius 3.
                // Condition (2) is more generous and easier to compute.
                const S maxSlope = 3 * sec;
                m0 = constrainSlope(m0, maxSlope);
                m1 = constrainSlope(m1, maxSlope);

                m0 = constrainSlope(m0, 3 * sec0);
                m1 = constrainSlope(m1, 3 * sec1);
            }

            const S t = (x - low->first) / interval;
            const S t2 = t * t;
            if (mCache) {
                // convert to cubic polynomial coefficients and compute
                m0 *= interval;
                m1 *= interval;
                const T dy = high->second - low->second;
                const S c0 = low->second;
                const S c1 = m0;
                const S c2 = 3 * dy - 2 * m0 - m1;
                const S c3 = m0 + m1 - 2 * dy;
                mMemo[low->first] = std::make_tuple(c1, c2, c3);
                return c0 + c1 * t + (c2 + c3 * t) * t2;
            } else {
                // classic Hermite interpolation
                const S t3 = t2 * t;
                const S h00 =  2 * t3 - 3 * t2     + 1;
                const S h10 =      t3 - 2 * t2 + t    ;
                const S h01 = -2 * t3 + 3 * t2        ;
                const S h11 =      t3     - t2        ;
                return h00 * low->second + (h10 * m0 + h11 * m1) * interval + h01 * high->second;
            }
        } // default
        }
    }

    InterpolatorType getInterpolatorType() const {
        return mInterpolatorType;
    }

    status_t setInterpolatorType(InterpolatorType interpolatorType) {
        switch (interpolatorType) {
        case INTERPOLATOR_TYPE_STEP:   // Not continuous
        case INTERPOLATOR_TYPE_LINEAR: // C0
        case INTERPOLATOR_TYPE_CUBIC:  // C1
        case INTERPOLATOR_TYPE_CUBIC_MONOTONIC: // C1 + other constraints
        // case INTERPOLATOR_TYPE_CUBIC_C2:
            mInterpolatorType = interpolatorType;
            return NO_ERROR;
        default:
            ALOGE("invalid interpolatorType: %d", interpolatorType);
            return BAD_VALUE;
        }
    }

    T getFirstSlope() const {
        return mFirstSlope;
    }

    void setFirstSlope(T slope) {
        mFirstSlope = slope;
    }

    T getLastSlope() const {
        return mLastSlope;
    }

    void setLastSlope(T slope) {
        mLastSlope = slope;
    }

    void clearCache() {
        mMemo.clear();
    }

    status_t writeToParcel(Parcel *parcel) const {
        if (parcel == nullptr) {
            return BAD_VALUE;
        }
        status_t res = parcel->writeInt32(mInterpolatorType)
                ?: parcel->writeFloat(mFirstSlope)
                ?: parcel->writeFloat(mLastSlope)
                ?: parcel->writeUint32((uint32_t)this->size()); // silent truncation
        if (res != NO_ERROR) {
            return res;
        }
        for (const auto &pt : *this) {
            res = parcel->writeFloat(pt.first)
                    ?: parcel->writeFloat(pt.second);
            if (res != NO_ERROR) {
                return res;
            }
        }
        return NO_ERROR;
    }

    status_t readFromParcel(const Parcel &parcel) {
        this->clear();
        int32_t type;
        uint32_t size;
        status_t res = parcel.readInt32(&type)
                        ?: parcel.readFloat(&mFirstSlope)
                        ?: parcel.readFloat(&mLastSlope)
                        ?: parcel.readUint32(&size)
                        ?: setInterpolatorType((InterpolatorType)type);
        if (res != NO_ERROR) {
            return res;
        }
        // Note: We don't need to check size is within some bounds as
        // the Parcel read will fail if size is incorrectly specified too large.
        float lastx;
        for (uint32_t i = 0; i < size; ++i) {
            float x, y;
            res = parcel.readFloat(&x)
                    ?: parcel.readFloat(&y);
            if (res != NO_ERROR) {
                return res;
            }
            if (i > 0 && !(x > lastx) /* handle nan */
                    || y != y /* handle nan */) {
                // This is a std::map object which imposes sorted order
                // automatically on emplace.
                // Nevertheless for reading from a Parcel,
                // we require that the points be specified monotonic in x.
                return BAD_VALUE;
            }
            this->emplace(x, y);
            lastx = x;
        }
        return NO_ERROR;
    }

    std::string toString() const {
        std::stringstream ss;
        ss << "Interpolator{mInterpolatorType=" << static_cast<int32_t>(mInterpolatorType);
        ss << ", mFirstSlope=" << mFirstSlope;
        ss << ", mLastSlope=" << mLastSlope;
        ss << ", {";
        bool first = true;
        for (const auto &pt : *this) {
            if (first) {
                first = false;
                ss << "{";
            } else {
                ss << ", {";
            }
            ss << pt.first << ", " << pt.second << "}";
        }
        ss << "}}";
        return ss.str();
    }

private:
    static S constrainSlope(S slope, S maxSlope) {
        if (maxSlope > 0) {
            slope = std::min(slope, maxSlope);
            slope = std::max(slope, S(0)); // not globally monotonic
        } else {
            slope = std::max(slope, maxSlope);
            slope = std::min(slope, S(0)); // not globally monotonic
        }
        return slope;
    }

    InterpolatorType mInterpolatorType;
    bool mCache; // whether we cache spline coefficient computation

    // for cubic interpolation, the boundary conditions in slope.
    S mFirstSlope;
    S mLastSlope;

    // spline cubic polynomial coefficient cache
    std::unordered_map<S, std::tuple<S /* c1 */, S /* c2 */, S /* c3 */>> mMemo;
}; // Interpolator

} // namespace android

#pragma pop_macro("LOG_TAG")

#endif // ANDROID_INTERPOLATOR_H