1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17package org.apache.commons.math.analysis.interpolation;
18
19import org.apache.commons.math.exception.DimensionMismatchException;
20import org.apache.commons.math.exception.NoDataException;
21import org.apache.commons.math.MathException;
22import org.apache.commons.math.util.MathUtils;
23
24/**
25 * Generates a tricubic interpolating function.
26 *
27 * @version $Revision$ $Date$
28 * @since 2.2
29 */
30public class TricubicSplineInterpolator
31    implements TrivariateRealGridInterpolator {
32    /**
33     * {@inheritDoc}
34     */
35    public TricubicSplineInterpolatingFunction interpolate(final double[] xval,
36                                                           final double[] yval,
37                                                           final double[] zval,
38                                                           final double[][][] fval)
39        throws MathException {
40        if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) {
41            throw new NoDataException();
42        }
43        if (xval.length != fval.length) {
44            throw new DimensionMismatchException(xval.length, fval.length);
45        }
46
47        MathUtils.checkOrder(xval);
48        MathUtils.checkOrder(yval);
49        MathUtils.checkOrder(zval);
50
51        final int xLen = xval.length;
52        final int yLen = yval.length;
53        final int zLen = zval.length;
54
55        // Samples, re-ordered as (z, x, y) and (y, z, x) tuplets
56        // fvalXY[k][i][j] = f(xval[i], yval[j], zval[k])
57        // fvalZX[j][k][i] = f(xval[i], yval[j], zval[k])
58        final double[][][] fvalXY = new double[zLen][xLen][yLen];
59        final double[][][] fvalZX = new double[yLen][zLen][xLen];
60        for (int i = 0; i < xLen; i++) {
61            if (fval[i].length != yLen) {
62                throw new DimensionMismatchException(fval[i].length, yLen);
63            }
64
65            for (int j = 0; j < yLen; j++) {
66                if (fval[i][j].length != zLen) {
67                    throw new DimensionMismatchException(fval[i][j].length, zLen);
68                }
69
70                for (int k = 0; k < zLen; k++) {
71                    final double v = fval[i][j][k];
72                    fvalXY[k][i][j] = v;
73                    fvalZX[j][k][i] = v;
74                }
75            }
76        }
77
78        final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator();
79
80        // For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z
81        final BicubicSplineInterpolatingFunction[] xSplineYZ
82            = new BicubicSplineInterpolatingFunction[xLen];
83        for (int i = 0; i < xLen; i++) {
84            xSplineYZ[i] = bsi.interpolate(yval, zval, fval[i]);
85        }
86
87        // For each line y[j] (0 <= j < yLen), construct a 2D spline in z and x
88        final BicubicSplineInterpolatingFunction[] ySplineZX
89            = new BicubicSplineInterpolatingFunction[yLen];
90        for (int j = 0; j < yLen; j++) {
91            ySplineZX[j] = bsi.interpolate(zval, xval, fvalZX[j]);
92        }
93
94        // For each line z[k] (0 <= k < zLen), construct a 2D spline in x and y
95        final BicubicSplineInterpolatingFunction[] zSplineXY
96            = new BicubicSplineInterpolatingFunction[zLen];
97        for (int k = 0; k < zLen; k++) {
98            zSplineXY[k] = bsi.interpolate(xval, yval, fvalXY[k]);
99        }
100
101        // Partial derivatives wrt x and wrt y
102        final double[][][] dFdX = new double[xLen][yLen][zLen];
103        final double[][][] dFdY = new double[xLen][yLen][zLen];
104        final double[][][] d2FdXdY = new double[xLen][yLen][zLen];
105        for (int k = 0; k < zLen; k++) {
106            final BicubicSplineInterpolatingFunction f = zSplineXY[k];
107            for (int i = 0; i < xLen; i++) {
108                final double x = xval[i];
109                for (int j = 0; j < yLen; j++) {
110                    final double y = yval[j];
111                    dFdX[i][j][k] = f.partialDerivativeX(x, y);
112                    dFdY[i][j][k] = f.partialDerivativeY(x, y);
113                    d2FdXdY[i][j][k] = f.partialDerivativeXY(x, y);
114                }
115            }
116        }
117
118        // Partial derivatives wrt y and wrt z
119        final double[][][] dFdZ = new double[xLen][yLen][zLen];
120        final double[][][] d2FdYdZ = new double[xLen][yLen][zLen];
121        for (int i = 0; i < xLen; i++) {
122            final BicubicSplineInterpolatingFunction f = xSplineYZ[i];
123            for (int j = 0; j < yLen; j++) {
124                final double y = yval[j];
125                for (int k = 0; k < zLen; k++) {
126                    final double z = zval[k];
127                    dFdZ[i][j][k] = f.partialDerivativeY(y, z);
128                    d2FdYdZ[i][j][k] = f.partialDerivativeXY(y, z);
129                }
130            }
131        }
132
133        // Partial derivatives wrt x and wrt z
134        final double[][][] d2FdZdX = new double[xLen][yLen][zLen];
135        for (int j = 0; j < yLen; j++) {
136            final BicubicSplineInterpolatingFunction f = ySplineZX[j];
137            for (int k = 0; k < zLen; k++) {
138                final double z = zval[k];
139                for (int i = 0; i < xLen; i++) {
140                    final double x = xval[i];
141                    d2FdZdX[i][j][k] = f.partialDerivativeXY(z, x);
142                }
143            }
144        }
145
146        // Third partial cross-derivatives
147        final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen];
148        for (int i = 0; i < xLen ; i++) {
149            final int nI = nextIndex(i, xLen);
150            final int pI = previousIndex(i);
151            for (int j = 0; j < yLen; j++) {
152                final int nJ = nextIndex(j, yLen);
153                final int pJ = previousIndex(j);
154                for (int k = 0; k < zLen; k++) {
155                    final int nK = nextIndex(k, zLen);
156                    final int pK = previousIndex(k);
157
158                    // XXX Not sure about this formula
159                    d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] -
160                                          fval[pI][nJ][nK] + fval[pI][pJ][nK] -
161                                          fval[nI][nJ][pK] + fval[nI][pJ][pK] +
162                                          fval[pI][nJ][pK] - fval[pI][pJ][pK]) /
163                        ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]) * (zval[nK] - zval[pK])) ;
164                }
165            }
166        }
167
168        // Create the interpolating splines
169        return new TricubicSplineInterpolatingFunction(xval, yval, zval, fval,
170                                                       dFdX, dFdY, dFdZ,
171                                                       d2FdXdY, d2FdZdX, d2FdYdZ,
172                                                       d3FdXdYdZ);
173    }
174
175    /**
176     * Compute the next index of an array, clipping if necessary.
177     * It is assumed (but not checked) that {@code i} is larger than or equal to 0}.
178     *
179     * @param i Index
180     * @param max Upper limit of the array
181     * @return the next index
182     */
183    private int nextIndex(int i, int max) {
184        final int index = i + 1;
185        return index < max ? index : index - 1;
186    }
187    /**
188     * Compute the previous index of an array, clipping if necessary.
189     * It is assumed (but not checked) that {@code i} is smaller than the size of the array.
190     *
191     * @param i Index
192     * @return the previous index
193     */
194    private int previousIndex(int i) {
195        final int index = i - 1;
196        return index >= 0 ? index : 0;
197    }
198}
199