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 }