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