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 3.4
029 */
030public class TricubicInterpolator
031    implements TrivariateGridInterpolator {
032    /**
033     * {@inheritDoc}
034     */
035    public TricubicInterpolatingFunction 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        // Approximation to the partial derivatives using finite differences.
057        final double[][][] dFdX = new double[xLen][yLen][zLen];
058        final double[][][] dFdY = new double[xLen][yLen][zLen];
059        final double[][][] dFdZ = new double[xLen][yLen][zLen];
060        final double[][][] d2FdXdY = new double[xLen][yLen][zLen];
061        final double[][][] d2FdXdZ = new double[xLen][yLen][zLen];
062        final double[][][] d2FdYdZ = new double[xLen][yLen][zLen];
063        final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen];
064
065        for (int i = 1; i < xLen - 1; i++) {
066            if (yval.length != fval[i].length) {
067                throw new DimensionMismatchException(yval.length, fval[i].length);
068            }
069
070            final int nI = i + 1;
071            final int pI = i - 1;
072
073            final double nX = xval[nI];
074            final double pX = xval[pI];
075
076            final double deltaX = nX - pX;
077
078            for (int j = 1; j < yLen - 1; j++) {
079                if (zval.length != fval[i][j].length) {
080                    throw new DimensionMismatchException(zval.length, fval[i][j].length);
081                }
082
083                final int nJ = j + 1;
084                final int pJ = j - 1;
085
086                final double nY = yval[nJ];
087                final double pY = yval[pJ];
088
089                final double deltaY = nY - pY;
090                final double deltaXY = deltaX * deltaY;
091
092                for (int k = 1; k < zLen - 1; k++) {
093                    final int nK = k + 1;
094                    final int pK = k - 1;
095
096                    final double nZ = zval[nK];
097                    final double pZ = zval[pK];
098
099                    final double deltaZ = nZ - pZ;
100
101                    dFdX[i][j][k] = (fval[nI][j][k] - fval[pI][j][k]) / deltaX;
102                    dFdY[i][j][k] = (fval[i][nJ][k] - fval[i][pJ][k]) / deltaY;
103                    dFdZ[i][j][k] = (fval[i][j][nK] - fval[i][j][pK]) / deltaZ;
104
105                    final double deltaXZ = deltaX * deltaZ;
106                    final double deltaYZ = deltaY * deltaZ;
107
108                    d2FdXdY[i][j][k] = (fval[nI][nJ][k] - fval[nI][pJ][k] - fval[pI][nJ][k] + fval[pI][pJ][k]) / deltaXY;
109                    d2FdXdZ[i][j][k] = (fval[nI][j][nK] - fval[nI][j][pK] - fval[pI][j][nK] + fval[pI][j][pK]) / deltaXZ;
110                    d2FdYdZ[i][j][k] = (fval[i][nJ][nK] - fval[i][nJ][pK] - fval[i][pJ][nK] + fval[i][pJ][pK]) / deltaYZ;
111
112                    final double deltaXYZ = deltaXY * deltaZ;
113
114                    d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] -
115                                          fval[pI][nJ][nK] + fval[pI][pJ][nK] -
116                                          fval[nI][nJ][pK] + fval[nI][pJ][pK] +
117                                          fval[pI][nJ][pK] - fval[pI][pJ][pK]) / deltaXYZ;
118                }
119            }
120        }
121
122        // Create the interpolating function.
123        return new TricubicInterpolatingFunction(xval, yval, zval, fval,
124                                                 dFdX, dFdY, dFdZ,
125                                                 d2FdXdY, d2FdXdZ, d2FdYdZ,
126                                                 d3FdXdYdZ) {
127            /** {@inheritDoc} */
128            @Override
129            public boolean isValidPoint(double x, double y, double z) {
130                if (x < xval[1] ||
131                    x > xval[xval.length - 2] ||
132                    y < yval[1] ||
133                    y > yval[yval.length - 2] ||
134                    z < zval[1] ||
135                    z > zval[zval.length - 2]) {
136                    return false;
137                } else {
138                    return true;
139                }
140            }
141        };
142    }
143}