    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.exception.DimensionMismatchException;
    020    import org.apache.commons.math3.exception.NoDataException;
    021    import org.apache.commons.math3.exception.NonMonotonicSequenceException;
    022    import org.apache.commons.math3.util.MathArrays;
    023
    024    /**
    025     * Generates a tricubic interpolating function.
    026     *
    027     * @since 2.2
    028     * @version $Id: TricubicSplineInterpolator.java 1379904 2012-09-01 23:54:52Z erans $
    029     */
    030    public class TricubicSplineInterpolator
    031        implements TrivariateGridInterpolator {
    032        /**
    033         * {@inheritDoc}
    034         */
    035        public TricubicSplineInterpolatingFunction interpolate(final double[] xval,
    036                                                               final double[] yval,
    037                                                               final double[] zval,
    038                                                               final double[][][] fval)
    039            throws NoDataException,
    040                   DimensionMismatchException,
    041                   NonMonotonicSequenceException {
    042            if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) {
    043                throw new NoDataException();
    044            }
    045            if (xval.length != fval.length) {
    046                throw new DimensionMismatchException(xval.length, fval.length);
    047            }
    048
    049            MathArrays.checkOrder(xval);
    050            MathArrays.checkOrder(yval);
    051            MathArrays.checkOrder(zval);
    052
    053            final int xLen = xval.length;
    054            final int yLen = yval.length;
    055            final int zLen = zval.length;
    056
    057            // Samples, re-ordered as (z, x, y) and (y, z, x) tuplets
    058            // fvalXY[k][i][j] = f(xval[i], yval[j], zval[k])
    059            // fvalZX[j][k][i] = f(xval[i], yval[j], zval[k])
    060            final double[][][] fvalXY = new double[zLen][xLen][yLen];
    061            final double[][][] fvalZX = new double[yLen][zLen][xLen];
    062            for (int i = 0; i < xLen; i++) {
    063                if (fval[i].length != yLen) {
    064                    throw new DimensionMismatchException(fval[i].length, yLen);
    065                }
    066
    067                for (int j = 0; j < yLen; j++) {
    068                    if (fval[i][j].length != zLen) {
    069                        throw new DimensionMismatchException(fval[i][j].length, zLen);
    070                    }
    071
    072                    for (int k = 0; k < zLen; k++) {
    073                        final double v = fval[i][j][k];
    074                        fvalXY[k][i][j] = v;
    075                        fvalZX[j][k][i] = v;
    076                    }
    077                }
    078            }
    079
    080            final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator();
    081
    082            // For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z
    083            final BicubicSplineInterpolatingFunction[] xSplineYZ
    084                = new BicubicSplineInterpolatingFunction[xLen];
    085            for (int i = 0; i < xLen; i++) {
    086                xSplineYZ[i] = bsi.interpolate(yval, zval, fval[i]);
    087            }
    088
    089            // For each line y[j] (0 <= j < yLen), construct a 2D spline in z and x
    090            final BicubicSplineInterpolatingFunction[] ySplineZX
    091                = new BicubicSplineInterpolatingFunction[yLen];
    092            for (int j = 0; j < yLen; j++) {
    093                ySplineZX[j] = bsi.interpolate(zval, xval, fvalZX[j]);
    094            }
    095
    096            // For each line z[k] (0 <= k < zLen), construct a 2D spline in x and y
    097            final BicubicSplineInterpolatingFunction[] zSplineXY
    098                = new BicubicSplineInterpolatingFunction[zLen];
    099            for (int k = 0; k < zLen; k++) {
    100                zSplineXY[k] = bsi.interpolate(xval, yval, fvalXY[k]);
    101            }
    102
    103            // Partial derivatives wrt x and wrt y
    104            final double[][][] dFdX = new double[xLen][yLen][zLen];
    105            final double[][][] dFdY = new double[xLen][yLen][zLen];
    106            final double[][][] d2FdXdY = new double[xLen][yLen][zLen];
    107            for (int k = 0; k < zLen; k++) {
    108                final BicubicSplineInterpolatingFunction f = zSplineXY[k];
    109                for (int i = 0; i < xLen; i++) {
    110                    final double x = xval[i];
    111                    for (int j = 0; j < yLen; j++) {
    112                        final double y = yval[j];
    113                        dFdX[i][j][k] = f.partialDerivativeX(x, y);
    114                        dFdY[i][j][k] = f.partialDerivativeY(x, y);
    115                        d2FdXdY[i][j][k] = f.partialDerivativeXY(x, y);
    116                    }
    117                }
    118            }
    119
    120            // Partial derivatives wrt y and wrt z
    121            final double[][][] dFdZ = new double[xLen][yLen][zLen];
    122            final double[][][] d2FdYdZ = new double[xLen][yLen][zLen];
    123            for (int i = 0; i < xLen; i++) {
    124                final BicubicSplineInterpolatingFunction f = xSplineYZ[i];
    125                for (int j = 0; j < yLen; j++) {
    126                    final double y = yval[j];
    127                    for (int k = 0; k < zLen; k++) {
    128                        final double z = zval[k];
    129                        dFdZ[i][j][k] = f.partialDerivativeY(y, z);
    130                        d2FdYdZ[i][j][k] = f.partialDerivativeXY(y, z);
    131                    }
    132                }
    133            }
    134
    135            // Partial derivatives wrt x and wrt z
    136            final double[][][] d2FdZdX = new double[xLen][yLen][zLen];
    137            for (int j = 0; j < yLen; j++) {
    138                final BicubicSplineInterpolatingFunction f = ySplineZX[j];
    139                for (int k = 0; k < zLen; k++) {
    140                    final double z = zval[k];
    141                    for (int i = 0; i < xLen; i++) {
    142                        final double x = xval[i];
    143                        d2FdZdX[i][j][k] = f.partialDerivativeXY(z, x);
    144                    }
    145                }
    146            }
    147
    148            // Third partial cross-derivatives
    149            final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen];
    150            for (int i = 0; i < xLen ; i++) {
    151                final int nI = nextIndex(i, xLen);
    152                final int pI = previousIndex(i);
    153                for (int j = 0; j < yLen; j++) {
    154                    final int nJ = nextIndex(j, yLen);
    155                    final int pJ = previousIndex(j);
    156                    for (int k = 0; k < zLen; k++) {
    157                        final int nK = nextIndex(k, zLen);
    158                        final int pK = previousIndex(k);
    159
    160                        // XXX Not sure about this formula
    161                        d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] -
    162                                              fval[pI][nJ][nK] + fval[pI][pJ][nK] -
    163                                              fval[nI][nJ][pK] + fval[nI][pJ][pK] +
    164                                              fval[pI][nJ][pK] - fval[pI][pJ][pK]) /
    165                            ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]) * (zval[nK] - zval[pK])) ;
    166                    }
    167                }
    168            }
    169
    170            // Create the interpolating splines
    171            return new TricubicSplineInterpolatingFunction(xval, yval, zval, fval,
    172                                                           dFdX, dFdY, dFdZ,
    173                                                           d2FdXdY, d2FdZdX, d2FdYdZ,
    174                                                           d3FdXdYdZ);
    175        }
    176
    177        /**
    178         * Compute the next index of an array, clipping if necessary.
    179         * It is assumed (but not checked) that {@code i} is larger than or equal to 0}.
    180         *
    181         * @param i Index
    182         * @param max Upper limit of the array
    183         * @return the next index
    184         */
    185        private int nextIndex(int i, int max) {
    186            final int index = i + 1;
    187            return index < max ? index : index - 1;
    188        }
    189        /**
    190         * Compute the previous index of an array, clipping if necessary.
    191         * It is assumed (but not checked) that {@code i} is smaller than the size of the array.
    192         *
    193         * @param i Index
    194         * @return the previous index
    195         */
    196        private int previousIndex(int i) {
    197            final int index = i - 1;
    198            return index >= 0 ? index : 0;
    199        }
    200    }