18921c28c7333ad2b4d34f013904ad4737044f366John Hofordpackage com.android.gallery3d.filtershow.filters;
28921c28c7333ad2b4d34f013904ad4737044f366John Hoford
38921c28c7333ad2b4d34f013904ad4737044f366John Hoford
48921c28c7333ad2b4d34f013904ad4737044f366John Hofordpublic class SplineMath {
58921c28c7333ad2b4d34f013904ad4737044f366John Hoford    double[][] mPoints = new double[6][2];
68921c28c7333ad2b4d34f013904ad4737044f366John Hoford    double[] mDerivatives;
78921c28c7333ad2b4d34f013904ad4737044f366John Hoford    SplineMath(int n) {
88921c28c7333ad2b4d34f013904ad4737044f366John Hoford        mPoints = new double[n][2];
98921c28c7333ad2b4d34f013904ad4737044f366John Hoford    }
108921c28c7333ad2b4d34f013904ad4737044f366John Hoford
118921c28c7333ad2b4d34f013904ad4737044f366John Hoford    public void setPoint(int index, double x, double y) {
128921c28c7333ad2b4d34f013904ad4737044f366John Hoford        mPoints[index][0] = x;
138921c28c7333ad2b4d34f013904ad4737044f366John Hoford        mPoints[index][1] = y;
148921c28c7333ad2b4d34f013904ad4737044f366John Hoford        mDerivatives = null;
158921c28c7333ad2b4d34f013904ad4737044f366John Hoford    }
168921c28c7333ad2b4d34f013904ad4737044f366John Hoford
178921c28c7333ad2b4d34f013904ad4737044f366John Hoford    public float[][] calculatetCurve(int n) {
188921c28c7333ad2b4d34f013904ad4737044f366John Hoford        float[][] curve = new float[n][2];
198921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double[][] points = new double[mPoints.length][2];
208921c28c7333ad2b4d34f013904ad4737044f366John Hoford        for (int i = 0; i < mPoints.length; i++) {
218921c28c7333ad2b4d34f013904ad4737044f366John Hoford
228921c28c7333ad2b4d34f013904ad4737044f366John Hoford            points[i][0] = mPoints[i][0];
238921c28c7333ad2b4d34f013904ad4737044f366John Hoford            points[i][1] = mPoints[i][1];
248921c28c7333ad2b4d34f013904ad4737044f366John Hoford
258921c28c7333ad2b4d34f013904ad4737044f366John Hoford        }
268921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double[] derivatives = solveSystem(points);
278921c28c7333ad2b4d34f013904ad4737044f366John Hoford        float start = (float) points[0][0];
288921c28c7333ad2b4d34f013904ad4737044f366John Hoford        float end = (float) (points[points.length - 1][0]);
298921c28c7333ad2b4d34f013904ad4737044f366John Hoford
308921c28c7333ad2b4d34f013904ad4737044f366John Hoford        curve[0][0] = (float) (points[0][0]);
318921c28c7333ad2b4d34f013904ad4737044f366John Hoford        curve[0][1] = (float) (points[0][1]);
328921c28c7333ad2b4d34f013904ad4737044f366John Hoford        int last = curve.length - 1;
338921c28c7333ad2b4d34f013904ad4737044f366John Hoford        curve[last][0] = (float) (points[points.length - 1][0]);
348921c28c7333ad2b4d34f013904ad4737044f366John Hoford        curve[last][1] = (float) (points[points.length - 1][1]);
358921c28c7333ad2b4d34f013904ad4737044f366John Hoford
368921c28c7333ad2b4d34f013904ad4737044f366John Hoford        for (int i = 0; i < curve.length; i++) {
378921c28c7333ad2b4d34f013904ad4737044f366John Hoford
388921c28c7333ad2b4d34f013904ad4737044f366John Hoford            double[] cur = null;
398921c28c7333ad2b4d34f013904ad4737044f366John Hoford            double[] next = null;
408921c28c7333ad2b4d34f013904ad4737044f366John Hoford            double x = start + i * (end - start) / (curve.length - 1);
418921c28c7333ad2b4d34f013904ad4737044f366John Hoford            int pivot = 0;
428921c28c7333ad2b4d34f013904ad4737044f366John Hoford            for (int j = 0; j < points.length - 1; j++) {
438921c28c7333ad2b4d34f013904ad4737044f366John Hoford                if (x >= points[j][0] && x <= points[j + 1][0]) {
448921c28c7333ad2b4d34f013904ad4737044f366John Hoford                    pivot = j;
458921c28c7333ad2b4d34f013904ad4737044f366John Hoford                }
468921c28c7333ad2b4d34f013904ad4737044f366John Hoford            }
478921c28c7333ad2b4d34f013904ad4737044f366John Hoford            cur = points[pivot];
488921c28c7333ad2b4d34f013904ad4737044f366John Hoford            next = points[pivot + 1];
498921c28c7333ad2b4d34f013904ad4737044f366John Hoford            if (x <= next[0]) {
508921c28c7333ad2b4d34f013904ad4737044f366John Hoford                double x1 = cur[0];
518921c28c7333ad2b4d34f013904ad4737044f366John Hoford                double x2 = next[0];
528921c28c7333ad2b4d34f013904ad4737044f366John Hoford                double y1 = cur[1];
538921c28c7333ad2b4d34f013904ad4737044f366John Hoford                double y2 = next[1];
548921c28c7333ad2b4d34f013904ad4737044f366John Hoford
558921c28c7333ad2b4d34f013904ad4737044f366John Hoford                // Use the second derivatives to apply the cubic spline
568921c28c7333ad2b4d34f013904ad4737044f366John Hoford                // equation:
578921c28c7333ad2b4d34f013904ad4737044f366John Hoford                double delta = (x2 - x1);
588921c28c7333ad2b4d34f013904ad4737044f366John Hoford                double delta2 = delta * delta;
598921c28c7333ad2b4d34f013904ad4737044f366John Hoford                double b = (x - x1) / delta;
608921c28c7333ad2b4d34f013904ad4737044f366John Hoford                double a = 1 - b;
618921c28c7333ad2b4d34f013904ad4737044f366John Hoford                double ta = a * y1;
628921c28c7333ad2b4d34f013904ad4737044f366John Hoford                double tb = b * y2;
638921c28c7333ad2b4d34f013904ad4737044f366John Hoford                double tc = (a * a * a - a) * derivatives[pivot];
648921c28c7333ad2b4d34f013904ad4737044f366John Hoford                double td = (b * b * b - b) * derivatives[pivot + 1];
658921c28c7333ad2b4d34f013904ad4737044f366John Hoford                double y = ta + tb + (delta2 / 6) * (tc + td);
668921c28c7333ad2b4d34f013904ad4737044f366John Hoford
678921c28c7333ad2b4d34f013904ad4737044f366John Hoford                curve[i][0] = (float) (x);
688921c28c7333ad2b4d34f013904ad4737044f366John Hoford                curve[i][1] = (float) (y);
698921c28c7333ad2b4d34f013904ad4737044f366John Hoford            } else {
708921c28c7333ad2b4d34f013904ad4737044f366John Hoford                curve[i][0] = (float) (next[0]);
718921c28c7333ad2b4d34f013904ad4737044f366John Hoford                curve[i][1] = (float) (next[1]);
728921c28c7333ad2b4d34f013904ad4737044f366John Hoford            }
738921c28c7333ad2b4d34f013904ad4737044f366John Hoford        }
748921c28c7333ad2b4d34f013904ad4737044f366John Hoford        return curve;
758921c28c7333ad2b4d34f013904ad4737044f366John Hoford    }
768921c28c7333ad2b4d34f013904ad4737044f366John Hoford
778921c28c7333ad2b4d34f013904ad4737044f366John Hoford    public double getValue(double x) {
788921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double[] cur = null;
798921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double[] next = null;
808921c28c7333ad2b4d34f013904ad4737044f366John Hoford        if (mDerivatives == null)
818921c28c7333ad2b4d34f013904ad4737044f366John Hoford            mDerivatives = solveSystem(mPoints);
828921c28c7333ad2b4d34f013904ad4737044f366John Hoford        int pivot = 0;
838921c28c7333ad2b4d34f013904ad4737044f366John Hoford        for (int j = 0; j < mPoints.length - 1; j++) {
848921c28c7333ad2b4d34f013904ad4737044f366John Hoford            pivot = j;
858921c28c7333ad2b4d34f013904ad4737044f366John Hoford            if (x <= mPoints[j][0]) {
868921c28c7333ad2b4d34f013904ad4737044f366John Hoford                break;
878921c28c7333ad2b4d34f013904ad4737044f366John Hoford            }
888921c28c7333ad2b4d34f013904ad4737044f366John Hoford        }
898921c28c7333ad2b4d34f013904ad4737044f366John Hoford        cur = mPoints[pivot];
908921c28c7333ad2b4d34f013904ad4737044f366John Hoford        next = mPoints[pivot + 1];
918921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double x1 = cur[0];
928921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double x2 = next[0];
938921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double y1 = cur[1];
948921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double y2 = next[1];
958921c28c7333ad2b4d34f013904ad4737044f366John Hoford
968921c28c7333ad2b4d34f013904ad4737044f366John Hoford        // Use the second derivatives to apply the cubic spline
978921c28c7333ad2b4d34f013904ad4737044f366John Hoford        // equation:
988921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double delta = (x2 - x1);
998921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double delta2 = delta * delta;
1008921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double b = (x - x1) / delta;
1018921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double a = 1 - b;
1028921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double ta = a * y1;
1038921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double tb = b * y2;
1048921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double tc = (a * a * a - a) * mDerivatives[pivot];
1058921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double td = (b * b * b - b) * mDerivatives[pivot + 1];
1068921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double y = ta + tb + (delta2 / 6) * (tc + td);
1078921c28c7333ad2b4d34f013904ad4737044f366John Hoford
1088921c28c7333ad2b4d34f013904ad4737044f366John Hoford        return y;
1098921c28c7333ad2b4d34f013904ad4737044f366John Hoford
1108921c28c7333ad2b4d34f013904ad4737044f366John Hoford    }
1118921c28c7333ad2b4d34f013904ad4737044f366John Hoford
1128921c28c7333ad2b4d34f013904ad4737044f366John Hoford    double[] solveSystem(double[][] points) {
1138921c28c7333ad2b4d34f013904ad4737044f366John Hoford        int n = points.length;
1148921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double[][] system = new double[n][3];
1158921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double[] result = new double[n]; // d
1168921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double[] solution = new double[n]; // returned coefficients
1178921c28c7333ad2b4d34f013904ad4737044f366John Hoford        system[0][1] = 1;
1188921c28c7333ad2b4d34f013904ad4737044f366John Hoford        system[n - 1][1] = 1;
1198921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double d6 = 1.0 / 6.0;
1208921c28c7333ad2b4d34f013904ad4737044f366John Hoford        double d3 = 1.0 / 3.0;
1218921c28c7333ad2b4d34f013904ad4737044f366John Hoford
1228921c28c7333ad2b4d34f013904ad4737044f366John Hoford        // let's create a tridiagonal matrix representing the
1238921c28c7333ad2b4d34f013904ad4737044f366John Hoford        // system, and apply the TDMA algorithm to solve it
1248921c28c7333ad2b4d34f013904ad4737044f366John Hoford        // (see http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm)
1258921c28c7333ad2b4d34f013904ad4737044f366John Hoford        for (int i = 1; i < n - 1; i++) {
1268921c28c7333ad2b4d34f013904ad4737044f366John Hoford            double deltaPrevX = points[i][0] - points[i - 1][0];
1278921c28c7333ad2b4d34f013904ad4737044f366John Hoford            double deltaX = points[i + 1][0] - points[i - 1][0];
1288921c28c7333ad2b4d34f013904ad4737044f366John Hoford            double deltaNextX = points[i + 1][0] - points[i][0];
1298921c28c7333ad2b4d34f013904ad4737044f366John Hoford            double deltaNextY = points[i + 1][1] - points[i][1];
1308921c28c7333ad2b4d34f013904ad4737044f366John Hoford            double deltaPrevY = points[i][1] - points[i - 1][1];
1318921c28c7333ad2b4d34f013904ad4737044f366John Hoford            system[i][0] = d6 * deltaPrevX; // a_i
1328921c28c7333ad2b4d34f013904ad4737044f366John Hoford            system[i][1] = d3 * deltaX; // b_i
1338921c28c7333ad2b4d34f013904ad4737044f366John Hoford            system[i][2] = d6 * deltaNextX; // c_i
1348921c28c7333ad2b4d34f013904ad4737044f366John Hoford            result[i] = (deltaNextY / deltaNextX) - (deltaPrevY / deltaPrevX); // d_i
1358921c28c7333ad2b4d34f013904ad4737044f366John Hoford        }
1368921c28c7333ad2b4d34f013904ad4737044f366John Hoford
1378921c28c7333ad2b4d34f013904ad4737044f366John Hoford        // Forward sweep
1388921c28c7333ad2b4d34f013904ad4737044f366John Hoford        for (int i = 1; i < n; i++) {
1398921c28c7333ad2b4d34f013904ad4737044f366John Hoford            // m = a_i/b_i-1
1408921c28c7333ad2b4d34f013904ad4737044f366John Hoford            double m = system[i][0] / system[i - 1][1];
1418921c28c7333ad2b4d34f013904ad4737044f366John Hoford            // b_i = b_i - m(c_i-1)
1428921c28c7333ad2b4d34f013904ad4737044f366John Hoford            system[i][1] = system[i][1] - m * system[i - 1][2];
1438921c28c7333ad2b4d34f013904ad4737044f366John Hoford            // d_i = d_i - m(d_i-1)
1448921c28c7333ad2b4d34f013904ad4737044f366John Hoford            result[i] = result[i] - m * result[i - 1];
1458921c28c7333ad2b4d34f013904ad4737044f366John Hoford        }
1468921c28c7333ad2b4d34f013904ad4737044f366John Hoford
1478921c28c7333ad2b4d34f013904ad4737044f366John Hoford        // Back substitution
1488921c28c7333ad2b4d34f013904ad4737044f366John Hoford        solution[n - 1] = result[n - 1] / system[n - 1][1];
1498921c28c7333ad2b4d34f013904ad4737044f366John Hoford        for (int i = n - 2; i >= 0; --i) {
1508921c28c7333ad2b4d34f013904ad4737044f366John Hoford            solution[i] = (result[i] - system[i][2] * solution[i + 1]) / system[i][1];
1518921c28c7333ad2b4d34f013904ad4737044f366John Hoford        }
1528921c28c7333ad2b4d34f013904ad4737044f366John Hoford        return solution;
1538921c28c7333ad2b4d34f013904ad4737044f366John Hoford    }
1548921c28c7333ad2b4d34f013904ad4737044f366John Hoford
1558921c28c7333ad2b4d34f013904ad4737044f366John Hoford    public static void main(String[] args) {
1568921c28c7333ad2b4d34f013904ad4737044f366John Hoford        SplineMath s = new SplineMath(10);
1578921c28c7333ad2b4d34f013904ad4737044f366John Hoford        for (int i = 0; i < 10; i++) {
1588921c28c7333ad2b4d34f013904ad4737044f366John Hoford            s.setPoint(i, i, i);
1598921c28c7333ad2b4d34f013904ad4737044f366John Hoford        }
1608921c28c7333ad2b4d34f013904ad4737044f366John Hoford        float[][] curve = s.calculatetCurve(40);
1618921c28c7333ad2b4d34f013904ad4737044f366John Hoford
1628921c28c7333ad2b4d34f013904ad4737044f366John Hoford        for (int j = 0; j < curve.length; j++) {
1638921c28c7333ad2b4d34f013904ad4737044f366John Hoford            System.out.println(curve[j][0] + "," + curve[j][1]);
1648921c28c7333ad2b4d34f013904ad4737044f366John Hoford        }
1658921c28c7333ad2b4d34f013904ad4737044f366John Hoford    }
1668921c28c7333ad2b4d34f013904ad4737044f366John Hoford}
167