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.math.analysis.interpolation;
018    
019    import org.apache.commons.math.exception.DimensionMismatchException;
020    import org.apache.commons.math.exception.NoDataException;
021    import org.apache.commons.math.MathException;
022    import org.apache.commons.math.util.MathUtils;
023    
024    /**
025     * Generates a tricubic interpolating function.
026     *
027     * @version $Revision$ $Date$
028     * @since 2.2
029     */
030    public class TricubicSplineInterpolator
031        implements TrivariateRealGridInterpolator {
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 MathException {
040            if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) {
041                throw new NoDataException();
042            }
043            if (xval.length != fval.length) {
044                throw new DimensionMismatchException(xval.length, fval.length);
045            }
046    
047            MathUtils.checkOrder(xval);
048            MathUtils.checkOrder(yval);
049            MathUtils.checkOrder(zval);
050    
051            final int xLen = xval.length;
052            final int yLen = yval.length;
053            final int zLen = zval.length;
054    
055            // Samples, re-ordered as (z, x, y) and (y, z, x) tuplets
056            // fvalXY[k][i][j] = f(xval[i], yval[j], zval[k])
057            // fvalZX[j][k][i] = f(xval[i], yval[j], zval[k])
058            final double[][][] fvalXY = new double[zLen][xLen][yLen];
059            final double[][][] fvalZX = new double[yLen][zLen][xLen];
060            for (int i = 0; i < xLen; i++) {
061                if (fval[i].length != yLen) {
062                    throw new DimensionMismatchException(fval[i].length, yLen);
063                }
064    
065                for (int j = 0; j < yLen; j++) {
066                    if (fval[i][j].length != zLen) {
067                        throw new DimensionMismatchException(fval[i][j].length, zLen);
068                    }
069    
070                    for (int k = 0; k < zLen; k++) {
071                        final double v = fval[i][j][k];
072                        fvalXY[k][i][j] = v;
073                        fvalZX[j][k][i] = v;
074                    }
075                }
076            }
077    
078            final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator();
079    
080            // For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z
081            final BicubicSplineInterpolatingFunction[] xSplineYZ
082                = new BicubicSplineInterpolatingFunction[xLen];
083            for (int i = 0; i < xLen; i++) {
084                xSplineYZ[i] = bsi.interpolate(yval, zval, fval[i]);
085            }
086    
087            // For each line y[j] (0 <= j < yLen), construct a 2D spline in z and x
088            final BicubicSplineInterpolatingFunction[] ySplineZX
089                = new BicubicSplineInterpolatingFunction[yLen];
090            for (int j = 0; j < yLen; j++) {
091                ySplineZX[j] = bsi.interpolate(zval, xval, fvalZX[j]);
092            }
093    
094            // For each line z[k] (0 <= k < zLen), construct a 2D spline in x and y
095            final BicubicSplineInterpolatingFunction[] zSplineXY
096                = new BicubicSplineInterpolatingFunction[zLen];
097            for (int k = 0; k < zLen; k++) {
098                zSplineXY[k] = bsi.interpolate(xval, yval, fvalXY[k]);
099            }
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    }