    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.math3.analysis.interpolation;
    018
    019    import org.apache.commons.math3.analysis.BivariateFunction;
    020    import org.apache.commons.math3.exception.DimensionMismatchException;
    021    import org.apache.commons.math3.exception.NoDataException;
    022    import org.apache.commons.math3.exception.OutOfRangeException;
    023    import org.apache.commons.math3.exception.NonMonotonicSequenceException;
    024    import org.apache.commons.math3.util.MathArrays;
    025
    026    /**
    027     * Function that implements the
    028     * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
    029     * bicubic spline interpolation</a>.
    030     *
    031     * @since 2.1
    032     * @version $Id: BicubicSplineInterpolatingFunction.java 1379904 2012-09-01 23:54:52Z erans $
    033     */
    034    public class BicubicSplineInterpolatingFunction
    035        implements BivariateFunction {
    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 BivariateFunction[][][] 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 NonMonotonicSequenceException if {@code x} or {@code y} are
    089         * 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                   NoDataException,
    100                   NonMonotonicSequenceException {
    101            final int xLen = x.length;
    102            final int yLen = y.length;
    103
    104            if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
    105                throw new NoDataException();
    106            }
    107            if (xLen != f.length) {
    108                throw new DimensionMismatchException(xLen, f.length);
    109            }
    110            if (xLen != dFdX.length) {
    111                throw new DimensionMismatchException(xLen, dFdX.length);
    112            }
    113            if (xLen != dFdY.length) {
    114                throw new DimensionMismatchException(xLen, dFdY.length);
    115            }
    116            if (xLen != d2FdXdY.length) {
    117                throw new DimensionMismatchException(xLen, d2FdXdY.length);
    118            }
    119
    120            MathArrays.checkOrder(x);
    121            MathArrays.checkOrder(y);
    122
    123            xval = x.clone();
    124            yval = y.clone();
    125
    126            final int lastI = xLen - 1;
    127            final int lastJ = yLen - 1;
    128            splines = new BicubicSplineFunction[lastI][lastJ];
    129
    130            for (int i = 0; i < lastI; i++) {
    131                if (f[i].length != yLen) {
    132                    throw new DimensionMismatchException(f[i].length, yLen);
    133                }
    134                if (dFdX[i].length != yLen) {
    135                    throw new DimensionMismatchException(dFdX[i].length, yLen);
    136                }
    137                if (dFdY[i].length != yLen) {
    138                    throw new DimensionMismatchException(dFdY[i].length, yLen);
    139                }
    140                if (d2FdXdY[i].length != yLen) {
    141                    throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
    142                }
    143                final int ip1 = i + 1;
    144                for (int j = 0; j < lastJ; j++) {
    145                    final int jp1 = j + 1;
    146                    final double[] beta = new double[] {
    147                        f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
    148                        dFdX[i][j], dFdX[ip1][j], dFdX[i][jp1], dFdX[ip1][jp1],
    149                        dFdY[i][j], dFdY[ip1][j], dFdY[i][jp1], dFdY[ip1][jp1],
    150                        d2FdXdY[i][j], d2FdXdY[ip1][j], d2FdXdY[i][jp1], d2FdXdY[ip1][jp1]
    151                    };
    152
    153                    splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta));
    154                }
    155            }
    156        }
    157
    158        /**
    159         * {@inheritDoc}
    160         */
    161        public double value(double x, double y)
    162            throws OutOfRangeException {
    163            final int i = searchIndex(x, xval);
    164            if (i == -1) {
    165                throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]);
    166            }
    167            final int j = searchIndex(y, yval);
    168            if (j == -1) {
    169                throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]);
    170            }
    171
    172            final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
    173            final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
    174
    175            return splines[i][j].value(xN, yN);
    176        }
    177
    178        /**
    179         * @param x x-coordinate.
    180         * @param y y-coordinate.
    181         * @return the value at point (x, y) of the first partial derivative with
    182         * respect to x.
    183         * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
    184         * the range defined by the boundary values of {@code xval} (resp.
    185         * {@code yval}).
    186         */
    187        public double partialDerivativeX(double x, double y)
    188            throws OutOfRangeException {
    189            return partialDerivative(0, x, y);
    190        }
    191        /**
    192         * @param x x-coordinate.
    193         * @param y y-coordinate.
    194         * @return the value at point (x, y) of the first partial derivative with
    195         * respect to y.
    196         * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
    197         * the range defined by the boundary values of {@code xval} (resp.
    198         * {@code yval}).
    199         */
    200        public double partialDerivativeY(double x, double y)
    201            throws OutOfRangeException {
    202            return partialDerivative(1, x, y);
    203        }
    204        /**
    205         * @param x x-coordinate.
    206         * @param y y-coordinate.
    207         * @return the value at point (x, y) of the second partial derivative with
    208         * respect to x.
    209         * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
    210         * the range defined by the boundary values of {@code xval} (resp.
    211         * {@code yval}).
    212         */
    213        public double partialDerivativeXX(double x, double y)
    214            throws OutOfRangeException {
    215            return partialDerivative(2, x, y);
    216        }
    217        /**
    218         * @param x x-coordinate.
    219         * @param y y-coordinate.
    220         * @return the value at point (x, y) of the second partial derivative with
    221         * respect to y.
    222         * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
    223         * the range defined by the boundary values of {@code xval} (resp.
    224         * {@code yval}).
    225         */
    226        public double partialDerivativeYY(double x, double y)
    227            throws OutOfRangeException {
    228            return partialDerivative(3, x, y);
    229        }
    230        /**
    231         * @param x x-coordinate.
    232         * @param y y-coordinate.
    233         * @return the value at point (x, y) of the second partial cross-derivative.
    234         * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
    235         * the range defined by the boundary values of {@code xval} (resp.
    236         * {@code yval}).
    237         */
    238        public double partialDerivativeXY(double x, double y)
    239            throws OutOfRangeException {
    240            return partialDerivative(4, x, y);
    241        }
    242
    243        /**
    244         * @param which First index in {@link #partialDerivatives}.
    245         * @param x x-coordinate.
    246         * @param y y-coordinate.
    247         * @return the value at point (x, y) of the selected partial derivative.
    248         * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
    249         * the range defined by the boundary values of {@code xval} (resp.
    250         * {@code yval}).
    251         */
    252        private double partialDerivative(int which, double x, double y)
    253            throws OutOfRangeException {
    254            if (partialDerivatives == null) {
    255                computePartialDerivatives();
    256            }
    257
    258            final int i = searchIndex(x, xval);
    259            if (i == -1) {
    260                throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]);
    261            }
    262            final int j = searchIndex(y, yval);
    263            if (j == -1) {
    264                throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]);
    265            }
    266
    267            final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
    268            final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
    269
    270            return partialDerivatives[which][i][j].value(xN, yN);
    271        }
    272
    273        /**
    274         * Compute all partial derivatives.
    275         */
    276        private void computePartialDerivatives() {
    277            final int lastI = xval.length - 1;
    278            final int lastJ = yval.length - 1;
    279            partialDerivatives = new BivariateFunction[5][lastI][lastJ];
    280
    281            for (int i = 0; i < lastI; i++) {
    282                for (int j = 0; j < lastJ; j++) {
    283                    final BicubicSplineFunction f = splines[i][j];
    284                    partialDerivatives[0][i][j] = f.partialDerivativeX();
    285                    partialDerivatives[1][i][j] = f.partialDerivativeY();
    286                    partialDerivatives[2][i][j] = f.partialDerivativeXX();
    287                    partialDerivatives[3][i][j] = f.partialDerivativeYY();
    288                    partialDerivatives[4][i][j] = f.partialDerivativeXY();
    289                }
    290            }
    291        }
    292
    293        /**
    294         * @param c Coordinate.
    295         * @param val Coordinate samples.
    296         * @return the index in {@code val} corresponding to the interval
    297         * containing {@code c}, or {@code -1} if {@code c} is out of the
    298         * range defined by the boundary values of {@code val}.
    299         */
    300        private int searchIndex(double c, double[] val) {
    301            if (c < val[0]) {
    302                return -1;
    303            }
    304
    305            final int max = val.length;
    306            for (int i = 1; i < max; i++) {
    307                if (c <= val[i]) {
    308                    return i - 1;
    309                }
    310            }
    311
    312            return -1;
    313        }
    314
    315        /**
    316         * Compute the spline coefficients from the list of function values and
    317         * function partial derivatives values at the four corners of a grid
    318         * element. They must be specified in the following order:
    319         * <ul>
    320         *  <li>f(0,0)</li>
    321         *  <li>f(1,0)</li>
    322         *  <li>f(0,1)</li>
    323         *  <li>f(1,1)</li>
    324         *  <li>f<sub>x</sub>(0,0)</li>
    325         *  <li>f<sub>x</sub>(1,0)</li>
    326         *  <li>f<sub>x</sub>(0,1)</li>
    327         *  <li>f<sub>x</sub>(1,1)</li>
    328         *  <li>f<sub>y</sub>(0,0)</li>
    329         *  <li>f<sub>y</sub>(1,0)</li>
    330         *  <li>f<sub>y</sub>(0,1)</li>
    331         *  <li>f<sub>y</sub>(1,1)</li>
    332         *  <li>f<sub>xy</sub>(0,0)</li>
    333         *  <li>f<sub>xy</sub>(1,0)</li>
    334         *  <li>f<sub>xy</sub>(0,1)</li>
    335         *  <li>f<sub>xy</sub>(1,1)</li>
    336         * </ul>
    337         * where the subscripts indicate the partial derivative with respect to
    338         * the corresponding variable(s).
    339         *
    340         * @param beta List of function values and function partial derivatives
    341         * values.
    342         * @return the spline coefficients.
    343         */
    344        private double[] computeSplineCoefficients(double[] beta) {
    345            final double[] a = new double[16];
    346
    347            for (int i = 0; i < 16; i++) {
    348                double result = 0;
    349                final double[] row = AINV[i];
    350                for (int j = 0; j < 16; j++) {
    351                    result += row[j] * beta[j];
    352                }
    353                a[i] = result;
    354            }
    355
    356            return a;
    357        }
    358    }
    359
    360    /**
    361     * 2D-spline function.
    362     *
    363     * @version $Id: BicubicSplineInterpolatingFunction.java 1379904 2012-09-01 23:54:52Z erans $
    364     */
    365    class BicubicSplineFunction
    366        implements BivariateFunction {
    367
    368        /** Number of points. */
    369        private static final short N = 4;
    370
    371        /** Coefficients */
    372        private final double[][] a;
    373
    374        /** First partial derivative along x. */
    375        private BivariateFunction partialDerivativeX;
    376
    377        /** First partial derivative along y. */
    378        private BivariateFunction partialDerivativeY;
    379
    380        /** Second partial derivative along x. */
    381        private BivariateFunction partialDerivativeXX;
    382
    383        /** Second partial derivative along y. */
    384        private BivariateFunction partialDerivativeYY;
    385
    386        /** Second crossed partial derivative. */
    387        private BivariateFunction partialDerivativeXY;
    388
    389        /**
    390         * Simple constructor.
    391         * @param a Spline coefficients
    392         */
    393        public BicubicSplineFunction(double[] a) {
    394            this.a = new double[N][N];
    395            for (int i = 0; i < N; i++) {
    396                for (int j = 0; j < N; j++) {
    397                    this.a[i][j] = a[i + N * j];
    398                }
    399            }
    400        }
    401
    402        /**
    403         * {@inheritDoc}
    404         */
    405        public double value(double x, double y) {
    406            if (x < 0 || x > 1) {
    407                throw new OutOfRangeException(x, 0, 1);
    408            }
    409            if (y < 0 || y > 1) {
    410                throw new OutOfRangeException(y, 0, 1);
    411            }
    412
    413            final double x2 = x * x;
    414            final double x3 = x2 * x;
    415            final double[] pX = {1, x, x2, x3};
    416
    417            final double y2 = y * y;
    418            final double y3 = y2 * y;
    419            final double[] pY = {1, y, y2, y3};
    420
    421            return apply(pX, pY, a);
    422        }
    423
    424        /**
    425         * Compute the value of the bicubic polynomial.
    426         *
    427         * @param pX Powers of the x-coordinate.
    428         * @param pY Powers of the y-coordinate.
    429         * @param coeff Spline coefficients.
    430         * @return the interpolated value.
    431         */
    432        private double apply(double[] pX, double[] pY, double[][] coeff) {
    433            double result = 0;
    434            for (int i = 0; i < N; i++) {
    435                for (int j = 0; j < N; j++) {
    436                    result += coeff[i][j] * pX[i] * pY[j];
    437                }
    438            }
    439
    440            return result;
    441        }
    442
    443        /**
    444         * @return the partial derivative wrt {@code x}.
    445         */
    446        public BivariateFunction partialDerivativeX() {
    447            if (partialDerivativeX == null) {
    448                computePartialDerivatives();
    449            }
    450
    451            return partialDerivativeX;
    452        }
    453        /**
    454         * @return the partial derivative wrt {@code y}.
    455         */
    456        public BivariateFunction partialDerivativeY() {
    457            if (partialDerivativeY == null) {
    458                computePartialDerivatives();
    459            }
    460
    461            return partialDerivativeY;
    462        }
    463        /**
    464         * @return the second partial derivative wrt {@code x}.
    465         */
    466        public BivariateFunction partialDerivativeXX() {
    467            if (partialDerivativeXX == null) {
    468                computePartialDerivatives();
    469            }
    470
    471            return partialDerivativeXX;
    472        }
    473        /**
    474         * @return the second partial derivative wrt {@code y}.
    475         */
    476        public BivariateFunction partialDerivativeYY() {
    477            if (partialDerivativeYY == null) {
    478                computePartialDerivatives();
    479            }
    480
    481            return partialDerivativeYY;
    482        }
    483        /**
    484         * @return the second partial cross-derivative.
    485         */
    486        public BivariateFunction partialDerivativeXY() {
    487            if (partialDerivativeXY == null) {
    488                computePartialDerivatives();
    489            }
    490
    491            return partialDerivativeXY;
    492        }
    493
    494        /**
    495         * Compute all partial derivatives functions.
    496         */
    497        private void computePartialDerivatives() {
    498            final double[][] aX = new double[N][N];
    499            final double[][] aY = new double[N][N];
    500            final double[][] aXX = new double[N][N];
    501            final double[][] aYY = new double[N][N];
    502            final double[][] aXY = new double[N][N];
    503
    504            for (int i = 0; i < N; i++) {
    505                for (int j = 0; j < N; j++) {
    506                    final double c = a[i][j];
    507                    aX[i][j] = i * c;
    508                    aY[i][j] = j * c;
    509                    aXX[i][j] = (i - 1) * aX[i][j];
    510                    aYY[i][j] = (j - 1) * aY[i][j];
    511                    aXY[i][j] = j * aX[i][j];
    512                }
    513            }
    514
    515            partialDerivativeX = new BivariateFunction() {
    516                    public double value(double x, double y)  {
    517                        final double x2 = x * x;
    518                        final double[] pX = {0, 1, x, x2};
    519
    520                        final double y2 = y * y;
    521                        final double y3 = y2 * y;
    522                        final double[] pY = {1, y, y2, y3};
    523
    524                        return apply(pX, pY, aX);
    525                    }
    526                };
    527            partialDerivativeY = new BivariateFunction() {
    528                    public double value(double x, double y)  {
    529                        final double x2 = x * x;
    530                        final double x3 = x2 * x;
    531                        final double[] pX = {1, x, x2, x3};
    532
    533                        final double y2 = y * y;
    534                        final double[] pY = {0, 1, y, y2};
    535
    536                        return apply(pX, pY, aY);
    537                    }
    538                };
    539            partialDerivativeXX = new BivariateFunction() {
    540                    public double value(double x, double y)  {
    541                        final double[] pX = {0, 0, 1, x};
    542
    543                        final double y2 = y * y;
    544                        final double y3 = y2 * y;
    545                        final double[] pY = {1, y, y2, y3};
    546
    547                        return apply(pX, pY, aXX);
    548                    }
    549                };
    550            partialDerivativeYY = new BivariateFunction() {
    551                    public double value(double x, double y)  {
    552                        final double x2 = x * x;
    553                        final double x3 = x2 * x;
    554                        final double[] pX = {1, x, x2, x3};
    555
    556                        final double[] pY = {0, 0, 1, y};
    557
    558                        return apply(pX, pY, aYY);
    559                    }
    560                };
    561            partialDerivativeXY = new BivariateFunction() {
    562                    public double value(double x, double y)  {
    563                        final double x2 = x * x;
    564                        final double[] pX = {0, 1, x, x2};
    565
    566                        final double y2 = y * y;
    567                        final double[] pY = {0, 1, y, y2};
    568
    569                        return apply(pX, pY, aXY);
    570                    }
    571                };
    572        }
    573    }