    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math.analysis.interpolation;
    018
    019    import org.apache.commons.math.DimensionMismatchException;
    020    import org.apache.commons.math.FunctionEvaluationException;
    021    import org.apache.commons.math.analysis.BivariateRealFunction;
    022    import org.apache.commons.math.exception.NoDataException;
    023    import org.apache.commons.math.exception.OutOfRangeException;
    024    import org.apache.commons.math.util.MathUtils;
    025
    026    /**
    027     * Function that implements the
    028     * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
    029     * bicubic spline interpolation</a>.
    030     *
    031     * @version $Revision$ $Date$
    032     * @since 2.1
    033     */
    034    public class BicubicSplineInterpolatingFunction
    035        implements BivariateRealFunction {
    036        /**
    037         * Matrix to compute the spline coefficients from the function values
    038         * and function derivatives values
    039         */
    040        private static final double[][] AINV = {
    041            { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
    042            { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
    043            { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
    044            { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
    045            { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
    046            { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
    047            { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
    048            { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
    049            { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
    050            { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
    051            { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
    052            { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
    053            { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
    054            { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
    055            { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
    056            { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
    057        };
    058
    059        /** Samples x-coordinates */
    060        private final double[] xval;
    061        /** Samples y-coordinates */
    062        private final double[] yval;
    063        /** Set of cubic splines patching the whole data grid */
    064        private final BicubicSplineFunction[][] splines;
    065        /**
    066         * Partial derivatives
    067         * The value of the first index determines the kind of derivatives:
    068         * 0 = first partial derivatives wrt x
    069         * 1 = first partial derivatives wrt y
    070         * 2 = second partial derivatives wrt x
    071         * 3 = second partial derivatives wrt y
    072         * 4 = cross partial derivatives
    073         */
    074        private BivariateRealFunction[][][] partialDerivatives = null;
    075
    076        /**
    077         * @param x Sample values of the x-coordinate, in increasing order.
    078         * @param y Sample values of the y-coordinate, in increasing order.
    079         * @param f Values of the function on every grid point.
    080         * @param dFdX Values of the partial derivative of function with respect
    081         * to x on every grid point.
    082         * @param dFdY Values of the partial derivative of function with respect
    083         * to y on every grid point.
    084         * @param d2FdXdY Values of the cross partial derivative of function on
    085         * every grid point.
    086         * @throws DimensionMismatchException if the various arrays do not contain
    087         * the expected number of elements.
    088         * @throws org.apache.commons.math.exception.NonMonotonousSequenceException
    089         * if {@code x} or {@code y} are not strictly increasing.
    090         * @throws NoDataException if any of the arrays has zero length.
    091         */
    092        public BicubicSplineInterpolatingFunction(double[] x,
    093                                                  double[] y,
    094                                                  double[][] f,
    095                                                  double[][] dFdX,
    096                                                  double[][] dFdY,
    097                                                  double[][] d2FdXdY)
    098            throws DimensionMismatchException {
    099            final int xLen = x.length;
    100            final int yLen = y.length;
    101
    102            if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
    103                throw new NoDataException();
    104            }
    105            if (xLen != f.length) {
    106                throw new DimensionMismatchException(xLen, f.length);
    107            }
    108            if (xLen != dFdX.length) {
    109                throw new DimensionMismatchException(xLen, dFdX.length);
    110            }
    111            if (xLen != dFdY.length) {
    112                throw new DimensionMismatchException(xLen, dFdY.length);
    113            }
    114            if (xLen != d2FdXdY.length) {
    115                throw new DimensionMismatchException(xLen, d2FdXdY.length);
    116            }
    117
    118            MathUtils.checkOrder(x);
    119            MathUtils.checkOrder(y);
    120
    121            xval = x.clone();
    122            yval = y.clone();
    123
    124            final int lastI = xLen - 1;
    125            final int lastJ = yLen - 1;
    126            splines = new BicubicSplineFunction[lastI][lastJ];
    127
    128            for (int i = 0; i < lastI; i++) {
    129                if (f[i].length != yLen) {
    130                    throw new DimensionMismatchException(f[i].length, yLen);
    131                }
    132                if (dFdX[i].length != yLen) {
    133                    throw new DimensionMismatchException(dFdX[i].length, yLen);
    134                }
    135                if (dFdY[i].length != yLen) {
    136                    throw new DimensionMismatchException(dFdY[i].length, yLen);
    137                }
    138                if (d2FdXdY[i].length != yLen) {
    139                    throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
    140                }
    141                final int ip1 = i + 1;
    142                for (int j = 0; j < lastJ; j++) {
    143                    final int jp1 = j + 1;
    144                    final double[] beta = new double[] {
    145                        f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
    146                        dFdX[i][j], dFdX[ip1][j], dFdX[i][jp1], dFdX[ip1][jp1],
    147                        dFdY[i][j], dFdY[ip1][j], dFdY[i][jp1], dFdY[ip1][jp1],
    148                        d2FdXdY[i][j], d2FdXdY[ip1][j], d2FdXdY[i][jp1], d2FdXdY[ip1][jp1]
    149                    };
    150
    151                    splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta));
    152                }
    153            }
    154        }
    155
    156        /**
    157         * {@inheritDoc}
    158         */
    159        public double value(double x, double y) {
    160            final int i = searchIndex(x, xval);
    161            if (i == -1) {
    162                throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]);
    163            }
    164            final int j = searchIndex(y, yval);
    165            if (j == -1) {
    166                throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]);
    167            }
    168
    169            final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
    170            final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
    171
    172            return splines[i][j].value(xN, yN);
    173        }
    174
    175        /**
    176         * @param x x-coordinate.
    177         * @param y y-coordinate.
    178         * @return the value at point (x, y) of the first partial derivative with
    179         * respect to x.
    180         * @since 2.2
    181         */
    182        public double partialDerivativeX(double x, double y) {
    183            return partialDerivative(0, x, y);
    184        }
    185        /**
    186         * @param x x-coordinate.
    187         * @param y y-coordinate.
    188         * @return the value at point (x, y) of the first partial derivative with
    189         * respect to y.
    190         * @since 2.2
    191         */
    192        public double partialDerivativeY(double x, double y) {
    193            return partialDerivative(1, x, y);
    194        }
    195        /**
    196         * @param x x-coordinate.
    197         * @param y y-coordinate.
    198         * @return the value at point (x, y) of the second partial derivative with
    199         * respect to x.
    200         * @since 2.2
    201         */
    202        public double partialDerivativeXX(double x, double y) {
    203            return partialDerivative(2, x, y);
    204        }
    205        /**
    206         * @param x x-coordinate.
    207         * @param y y-coordinate.
    208         * @return the value at point (x, y) of the second partial derivative with
    209         * respect to y.
    210         * @since 2.2
    211         */
    212        public double partialDerivativeYY(double x, double y) {
    213            return partialDerivative(3, x, y);
    214        }
    215        /**
    216         * @param x x-coordinate.
    217         * @param y y-coordinate.
    218         * @return the value at point (x, y) of the second partial cross-derivative.
    219         * @since 2.2
    220         */
    221        public double partialDerivativeXY(double x, double y) {
    222            return partialDerivative(4, x, y);
    223        }
    224
    225        /**
    226         * @param which First index in {@link #partialDerivatives}.
    227         * @param x x-coordinate.
    228         * @param y y-coordinate.
    229         * @return the value at point (x, y) of the selected partial derivative.
    230         * @throws FunctionEvaluationException
    231         */
    232        private double partialDerivative(int which, double x, double y) {
    233            if (partialDerivatives == null) {
    234                computePartialDerivatives();
    235            }
    236
    237            final int i = searchIndex(x, xval);
    238            if (i == -1) {
    239                throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]);
    240            }
    241            final int j = searchIndex(y, yval);
    242            if (j == -1) {
    243                throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]);
    244            }
    245
    246            final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
    247            final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
    248
    249            try {
    250                return partialDerivatives[which][i][j].value(xN, yN);
    251            } catch (FunctionEvaluationException fee) {
    252                // this should never happen
    253                throw new RuntimeException(fee);
    254            }
    255
    256        }
    257
    258        /**
    259         * Compute all partial derivatives.
    260         */
    261        private void computePartialDerivatives() {
    262            final int lastI = xval.length - 1;
    263            final int lastJ = yval.length - 1;
    264            partialDerivatives = new BivariateRealFunction[5][lastI][lastJ];
    265
    266            for (int i = 0; i < lastI; i++) {
    267                for (int j = 0; j < lastJ; j++) {
    268                    final BicubicSplineFunction f = splines[i][j];
    269                    partialDerivatives[0][i][j] = f.partialDerivativeX();
    270                    partialDerivatives[1][i][j] = f.partialDerivativeY();
    271                    partialDerivatives[2][i][j] = f.partialDerivativeXX();
    272                    partialDerivatives[3][i][j] = f.partialDerivativeYY();
    273                    partialDerivatives[4][i][j] = f.partialDerivativeXY();
    274                }
    275            }
    276        }
    277
    278        /**
    279         * @param c Coordinate.
    280         * @param val Coordinate samples.
    281         * @return the index in {@code val} corresponding to the interval
    282         * containing {@code c}, or {@code -1} if {@code c} is out of the
    283         * range defined by the end values of {@code val}.
    284         */
    285        private int searchIndex(double c, double[] val) {
    286            if (c < val[0]) {
    287                return -1;
    288            }
    289
    290            final int max = val.length;
    291            for (int i = 1; i < max; i++) {
    292                if (c <= val[i]) {
    293                    return i - 1;
    294                }
    295            }
    296
    297            return -1;
    298        }
    299
    300        /**
    301         * Compute the spline coefficients from the list of function values and
    302         * function partial derivatives values at the four corners of a grid
    303         * element. They must be specified in the following order:
    304         * <ul>
    305         *  <li>f(0,0)</li>
    306         *  <li>f(1,0)</li>
    307         *  <li>f(0,1)</li>
    308         *  <li>f(1,1)</li>
    309         *  <li>f<sub>x</sub>(0,0)</li>
    310         *  <li>f<sub>x</sub>(1,0)</li>
    311         *  <li>f<sub>x</sub>(0,1)</li>
    312         *  <li>f<sub>x</sub>(1,1)</li>
    313         *  <li>f<sub>y</sub>(0,0)</li>
    314         *  <li>f<sub>y</sub>(1,0)</li>
    315         *  <li>f<sub>y</sub>(0,1)</li>
    316         *  <li>f<sub>y</sub>(1,1)</li>
    317         *  <li>f<sub>xy</sub>(0,0)</li>
    318         *  <li>f<sub>xy</sub>(1,0)</li>
    319         *  <li>f<sub>xy</sub>(0,1)</li>
    320         *  <li>f<sub>xy</sub>(1,1)</li>
    321         * </ul>
    322         * where the subscripts indicate the partial derivative with respect to
    323         * the corresponding variable(s).
    324         *
    325         * @param beta List of function values and function partial derivatives
    326         * values.
    327         * @return the spline coefficients.
    328         */
    329        private double[] computeSplineCoefficients(double[] beta) {
    330            final double[] a = new double[16];
    331
    332            for (int i = 0; i < 16; i++) {
    333                double result = 0;
    334                final double[] row = AINV[i];
    335                for (int j = 0; j < 16; j++) {
    336                    result += row[j] * beta[j];
    337                }
    338                a[i] = result;
    339            }
    340
    341            return a;
    342        }
    343    }
    344
    345    /**
    346     * 2D-spline function.
    347     *
    348     * @version $Revision$ $Date$
    349     */
    350    class BicubicSplineFunction
    351        implements BivariateRealFunction {
    352
    353        /** Number of points. */
    354        private static final short N = 4;
    355
    356        /** Coefficients */
    357        private final double[][] a;
    358
    359        /** First partial derivative along x. */
    360        private BivariateRealFunction partialDerivativeX;
    361
    362        /** First partial derivative along y. */
    363        private BivariateRealFunction partialDerivativeY;
    364
    365        /** Second partial derivative along x. */
    366        private BivariateRealFunction partialDerivativeXX;
    367
    368        /** Second partial derivative along y. */
    369        private BivariateRealFunction partialDerivativeYY;
    370
    371        /** Second crossed partial derivative. */
    372        private BivariateRealFunction partialDerivativeXY;
    373
    374        /**
    375         * Simple constructor.
    376         * @param a Spline coefficients
    377         */
    378        public BicubicSplineFunction(double[] a) {
    379            this.a = new double[N][N];
    380            for (int i = 0; i < N; i++) {
    381                for (int j = 0; j < N; j++) {
    382                    this.a[i][j] = a[i + N * j];
    383                }
    384            }
    385        }
    386
    387        /**
    388         * {@inheritDoc}
    389         */
    390        public double value(double x, double y) {
    391            if (x < 0 || x > 1) {
    392                throw new OutOfRangeException(x, 0, 1);
    393            }
    394            if (y < 0 || y > 1) {
    395                throw new OutOfRangeException(y, 0, 1);
    396            }
    397
    398            final double x2 = x * x;
    399            final double x3 = x2 * x;
    400            final double[] pX = {1, x, x2, x3};
    401
    402            final double y2 = y * y;
    403            final double y3 = y2 * y;
    404            final double[] pY = {1, y, y2, y3};
    405
    406            return apply(pX, pY, a);
    407        }
    408
    409        /**
    410         * Compute the value of the bicubic polynomial.
    411         *
    412         * @param pX Powers of the x-coordinate.
    413         * @param pY Powers of the y-coordinate.
    414         * @param coeff Spline coefficients.
    415         * @return the interpolated value.
    416         */
    417        private double apply(double[] pX, double[] pY, double[][] coeff) {
    418            double result = 0;
    419            for (int i = 0; i < N; i++) {
    420                for (int j = 0; j < N; j++) {
    421                    result += coeff[i][j] * pX[i] * pY[j];
    422                }
    423            }
    424
    425            return result;
    426        }
    427
    428        /**
    429         * @return the partial derivative wrt {@code x}.
    430         */
    431        public BivariateRealFunction partialDerivativeX() {
    432            if (partialDerivativeX == null) {
    433                computePartialDerivatives();
    434            }
    435
    436            return partialDerivativeX;
    437        }
    438        /**
    439         * @return the partial derivative wrt {@code y}.
    440         */
    441        public BivariateRealFunction partialDerivativeY() {
    442            if (partialDerivativeY == null) {
    443                computePartialDerivatives();
    444            }
    445
    446            return partialDerivativeY;
    447        }
    448        /**
    449         * @return the second partial derivative wrt {@code x}.
    450         */
    451        public BivariateRealFunction partialDerivativeXX() {
    452            if (partialDerivativeXX == null) {
    453                computePartialDerivatives();
    454            }
    455
    456            return partialDerivativeXX;
    457        }
    458        /**
    459         * @return the second partial derivative wrt {@code y}.
    460         */
    461        public BivariateRealFunction partialDerivativeYY() {
    462            if (partialDerivativeYY == null) {
    463                computePartialDerivatives();
    464            }
    465
    466            return partialDerivativeYY;
    467        }
    468        /**
    469         * @return the second partial cross-derivative.
    470         */
    471        public BivariateRealFunction partialDerivativeXY() {
    472            if (partialDerivativeXY == null) {
    473                computePartialDerivatives();
    474            }
    475
    476            return partialDerivativeXY;
    477        }
    478
    479        /**
    480         * Compute all partial derivatives functions.
    481         */
    482        private void computePartialDerivatives() {
    483            final double[][] aX = new double[N][N];
    484            final double[][] aY = new double[N][N];
    485            final double[][] aXX = new double[N][N];
    486            final double[][] aYY = new double[N][N];
    487            final double[][] aXY = new double[N][N];
    488
    489            for (int i = 0; i < N; i++) {
    490                for (int j = 0; j < N; j++) {
    491                    final double c = a[i][j];
    492                    aX[i][j] = i * c;
    493                    aY[i][j] = j * c;
    494                    aXX[i][j] = (i - 1) * aX[i][j];
    495                    aYY[i][j] = (j - 1) * aY[i][j];
    496                    aXY[i][j] = j * aX[i][j];
    497                }
    498            }
    499
    500            partialDerivativeX = new BivariateRealFunction() {
    501                    public double value(double x, double y)  {
    502                        final double x2 = x * x;
    503                        final double[] pX = {0, 1, x, x2};
    504
    505                        final double y2 = y * y;
    506                        final double y3 = y2 * y;
    507                        final double[] pY = {1, y, y2, y3};
    508
    509                        return apply(pX, pY, aX);
    510                    }
    511                };
    512            partialDerivativeY = new BivariateRealFunction() {
    513                    public double value(double x, double y)  {
    514                        final double x2 = x * x;
    515                        final double x3 = x2 * x;
    516                        final double[] pX = {1, x, x2, x3};
    517
    518                        final double y2 = y * y;
    519                        final double[] pY = {0, 1, y, y2};
    520
    521                        return apply(pX, pY, aY);
    522                    }
    523                };
    524            partialDerivativeXX = new BivariateRealFunction() {
    525                    public double value(double x, double y)  {
    526                        final double[] pX = {0, 0, 1, x};
    527
    528                        final double y2 = y * y;
    529                        final double y3 = y2 * y;
    530                        final double[] pY = {1, y, y2, y3};
    531
    532                        return apply(pX, pY, aXX);
    533                    }
    534                };
    535            partialDerivativeYY = new BivariateRealFunction() {
    536                    public double value(double x, double y)  {
    537                        final double x2 = x * x;
    538                        final double x3 = x2 * x;
    539                        final double[] pX = {1, x, x2, x3};
    540
    541                        final double[] pY = {0, 0, 1, y};
    542
    543                        return apply(pX, pY, aYY);
    544                    }
    545                };
    546            partialDerivativeXY = new BivariateRealFunction() {
    547                    public double value(double x, double y)  {
    548                        final double x2 = x * x;
    549                        final double[] pX = {0, 1, x, x2};
    550
    551                        final double y2 = y * y;
    552                        final double[] pY = {0, 1, y, y2};
    553
    554                        return apply(pX, pY, aXY);
    555                    }
    556                };
    557        }
    558    }