11a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown/*
21a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown * Copyright (C) 2012 The Android Open Source Project
31a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown *
41a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown * Licensed under the Apache License, Version 2.0 (the "License");
51a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown * you may not use this file except in compliance with the License.
61a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown * You may obtain a copy of the License at
71a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown *
81a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown *      http://www.apache.org/licenses/LICENSE-2.0
91a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown *
101a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown * Unless required by applicable law or agreed to in writing, software
111a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown * distributed under the License is distributed on an "AS IS" BASIS,
121a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
131a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown * See the License for the specific language governing permissions and
141a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown * limitations under the License.
151a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown */
161a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown
171a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brownpackage android.util;
181a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown
191a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown/**
201a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown * Performs spline interpolation given a set of control points.
211a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown * @hide
221a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown */
231a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brownpublic final class Spline {
241a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    private final float[] mX;
251a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    private final float[] mY;
261a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    private final float[] mM;
271a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown
281a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    private Spline(float[] x, float[] y, float[] m) {
291a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        mX = x;
301a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        mY = y;
311a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        mM = m;
321a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    }
331a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown
341a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    /**
351a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * Creates a monotone cubic spline from a given set of control points.
361a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     *
371a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * The spline is guaranteed to pass through each control point exactly.
381a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * Moreover, assuming the control points are monotonic (Y is non-decreasing or
391a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * non-increasing) then the interpolated values will also be monotonic.
401a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     *
411a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * This function uses the Fritsch-Carlson method for computing the spline parameters.
421a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * http://en.wikipedia.org/wiki/Monotone_cubic_interpolation
431a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     *
441a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * @param x The X component of the control points, strictly increasing.
451a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * @param y The Y component of the control points, monotonic.
461a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * @return
471a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     *
481a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * @throws IllegalArgumentException if the X or Y arrays are null, have
491a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * different lengths or have fewer than 2 values.
501a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * @throws IllegalArgumentException if the control points are not monotonic.
511a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     */
521a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    public static Spline createMonotoneCubicSpline(float[] x, float[] y) {
531a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        if (x == null || y == null || x.length != y.length || x.length < 2) {
541a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            throw new IllegalArgumentException("There must be at least two control "
551a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                    + "points and the arrays must be of equal length.");
561a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        }
571a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown
581a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        final int n = x.length;
591a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        float[] d = new float[n - 1]; // could optimize this out
601a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        float[] m = new float[n];
611a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown
621a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        // Compute slopes of secant lines between successive points.
631a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        for (int i = 0; i < n - 1; i++) {
641a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            float h = x[i + 1] - x[i];
651a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            if (h <= 0f) {
661a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                throw new IllegalArgumentException("The control points must all "
671a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                        + "have strictly increasing X values.");
681a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            }
691a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            d[i] = (y[i + 1] - y[i]) / h;
701a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        }
711a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown
721a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        // Initialize the tangents as the average of the secants.
731a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        m[0] = d[0];
741a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        for (int i = 1; i < n - 1; i++) {
751a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            m[i] = (d[i - 1] + d[i]) * 0.5f;
761a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        }
771a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        m[n - 1] = d[n - 2];
781a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown
791a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        // Update the tangents to preserve monotonicity.
801a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        for (int i = 0; i < n - 1; i++) {
811a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            if (d[i] == 0f) { // successive Y values are equal
821a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                m[i] = 0f;
831a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                m[i + 1] = 0f;
841a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            } else {
851a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                float a = m[i] / d[i];
861a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                float b = m[i + 1] / d[i];
871a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                if (a < 0f || b < 0f) {
881a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                    throw new IllegalArgumentException("The control points must have "
891a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                            + "monotonic Y values.");
901a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                }
911a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                float h = FloatMath.hypot(a, b);
921a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                if (h > 9f) {
931a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                    float t = 3f / h;
941a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                    m[i] = t * a * d[i];
951a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                    m[i + 1] = t * b * d[i];
961a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                }
971a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            }
981a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        }
991a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        return new Spline(x, y, m);
1001a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    }
1011a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown
1021a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    /**
1031a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * Interpolates the value of Y = f(X) for given X.
1041a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * Clamps X to the domain of the spline.
1051a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     *
1061a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * @param x The X value.
1071a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     * @return The interpolated Y = f(X) value.
1081a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown     */
1091a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    public float interpolate(float x) {
1101a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        // Handle the boundary cases.
1111a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        final int n = mX.length;
1121a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        if (Float.isNaN(x)) {
1131a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            return x;
1141a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        }
1151a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        if (x <= mX[0]) {
1161a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            return mY[0];
1171a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        }
1181a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        if (x >= mX[n - 1]) {
1191a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            return mY[n - 1];
1201a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        }
1211a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown
1221a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        // Find the index 'i' of the last point with smaller X.
1231a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        // We know this will be within the spline due to the boundary tests.
1241a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        int i = 0;
1251a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        while (x >= mX[i + 1]) {
1261a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            i += 1;
1271a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            if (x == mX[i]) {
1281a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                return mY[i];
1291a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            }
1301a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        }
1311a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown
1321a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        // Perform cubic Hermite spline interpolation.
1331a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        float h = mX[i + 1] - mX[i];
1341a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        float t = (x - mX[i]) / h;
1351a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        return (mY[i] * (1 + 2 * t) + h * mM[i] * t) * (1 - t) * (1 - t)
1361a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                + (mY[i + 1] * (3 - 2 * t) + h * mM[i + 1] * (t - 1)) * t * t;
1371a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    }
1381a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown
1391a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    // For debugging.
1401a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    @Override
1411a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    public String toString() {
1421a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        StringBuilder str = new StringBuilder();
1431a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        final int n = mX.length;
1441a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        str.append("[");
1451a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        for (int i = 0; i < n; i++) {
1461a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            if (i != 0) {
1471a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown                str.append(", ");
1481a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            }
1491a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            str.append("(").append(mX[i]);
1501a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            str.append(", ").append(mY[i]);
1511a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown            str.append(": ").append(mM[i]).append(")");
1521a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        }
1531a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        str.append("]");
1541a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown        return str.toString();
1551a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown    }
1561a30b55036c2279d72ba69cb1107ec5f6f40d5e9Jeff Brown}
157