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.util.MathUtils;
020 import org.apache.commons.math.MathRuntimeException;
021 import org.apache.commons.math.DimensionMismatchException;
022 import org.apache.commons.math.analysis.BivariateRealFunction;
023
024 /**
025 * Function that implements the
026 * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
027 * bicubic spline interpolation</a>.
028 *
029 * @version $Revision$ $Date$
030 * @since 2.1
031 */
032 public class BicubicSplineInterpolatingFunction
033 implements BivariateRealFunction {
034 /**
035 * Matrix to compute the spline coefficients from the function values
036 * and function derivatives values
037 */
038 private static final double[][] AINV = {
039 { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
040 { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
041 { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
042 { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
043 { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
044 { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
045 { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
046 { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
047 { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
048 { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
049 { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
050 { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
051 { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
052 { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
053 { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
054 { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
055 };
056
057 /** Samples x-coordinates */
058 private final double[] xval;
059 /** Samples y-coordinates */
060 private final double[] yval;
061 /** Set of cubic splines pacthing the whole data grid */
062 private final BicubicSplineFunction[][] splines;
063
064 /**
065 * @param x Sample values of the x-coordinate, in increasing order
066 * @param y Sample values of the y-coordinate, in increasing order
067 * @param z Values of the function on every grid point
068 * @param dZdX Values of the partial derivative of function with respect
069 * to x on every grid point
070 * @param dZdY Values of the partial derivative of function with respect
071 * to y on every grid point
072 * @param dZdXdY Values of the cross partial derivative of function on
073 * every grid point
074 * @throws DimensionMismatchException if the various arrays do not contain
075 * the expected number of elements.
076 * @throws IllegalArgumentException if {@code x} or {@code y} are not strictly
077 * increasing.
078 */
079 public BicubicSplineInterpolatingFunction(double[] x,
080 double[] y,
081 double[][] z,
082 double[][] dZdX,
083 double[][] dZdY,
084 double[][] dZdXdY)
085 throws DimensionMismatchException {
086 final int xLen = x.length;
087 final int yLen = y.length;
088
089 if (xLen == 0 || yLen == 0 || z.length == 0 || z[0].length == 0) {
090 throw MathRuntimeException.createIllegalArgumentException("no data");
091 }
092 if (xLen != z.length) {
093 throw new DimensionMismatchException(xLen, z.length);
094 }
095 if (xLen != dZdX.length) {
096 throw new DimensionMismatchException(xLen, dZdX.length);
097 }
098 if (xLen != dZdY.length) {
099 throw new DimensionMismatchException(xLen, dZdY.length);
100 }
101 if (xLen != dZdXdY.length) {
102 throw new DimensionMismatchException(xLen, dZdXdY.length);
103 }
104
105 MathUtils.checkOrder(x, 1, true);
106 MathUtils.checkOrder(y, 1, true);
107
108 xval = x.clone();
109 yval = y.clone();
110
111 final int lastI = xLen - 1;
112 final int lastJ = yLen - 1;
113 splines = new BicubicSplineFunction[lastI][lastJ];
114
115 for (int i = 0; i < lastI; i++) {
116 if (z[i].length != yLen) {
117 throw new DimensionMismatchException(z[i].length, yLen);
118 }
119 if (dZdX[i].length != yLen) {
120 throw new DimensionMismatchException(dZdX[i].length, yLen);
121 }
122 if (dZdY[i].length != yLen) {
123 throw new DimensionMismatchException(dZdY[i].length, yLen);
124 }
125 if (dZdXdY[i].length != yLen) {
126 throw new DimensionMismatchException(dZdXdY[i].length, yLen);
127 }
128 final int ip1 = i + 1;
129 for (int j = 0; j < lastJ; j++) {
130 final int jp1 = j + 1;
131 final double[] beta = new double[] {
132 z[i][j], z[ip1][j], z[i][jp1], z[ip1][jp1],
133 dZdX[i][j], dZdX[ip1][j], dZdX[i][jp1], dZdX[ip1][jp1],
134 dZdY[i][j], dZdY[ip1][j], dZdY[i][jp1], dZdY[ip1][jp1],
135 dZdXdY[i][j], dZdXdY[ip1][j], dZdXdY[i][jp1], dZdXdY[ip1][jp1]
136 };
137
138 splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta));
139 }
140 }
141 }
142
143 /**
144 * {@inheritDoc}
145 */
146 public double value(double x, double y) {
147 final int i = searchIndex(x, xval);
148 if (i == -1) {
149 throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
150 x, xval[0], xval[xval.length - 1]);
151 }
152 final int j = searchIndex(y, yval);
153 if (j == -1) {
154 throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
155 y, yval[0], yval[yval.length - 1]);
156 }
157
158 final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
159 final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
160
161 return splines[i][j].value(xN, yN);
162 }
163
164 /**
165 * @param c coordinate
166 * @param val coordinate samples
167 * @return the index in {@code val} corresponding to the interval
168 * containing {@code c}, or {@code -1} if {@code c} is out of the
169 * range defined by the end values of {@code val}
170 */
171 private int searchIndex(double c, double[] val) {
172 if (c < val[0]) {
173 return -1;
174 }
175
176 int max = val.length;
177 for (int i = 1; i < max; i++) {
178 if (c <= val[i]) {
179 return i - 1;
180 }
181 }
182
183 return -1;
184 }
185
186 /**
187 * Compute the spline coefficients from the list of function values and
188 * function partial derivatives values at the four corners of a grid
189 * element. They must be specified in the following order:
190 * <ul>
191 * <li>f(0,0)</li>
192 * <li>f(1,0)</li>
193 * <li>f(0,1)</li>
194 * <li>f(1,1)</li>
195 * <li>fx(0,0)</li>
196 * <li>fx(1,0)</li>
197 * <li>fx(0,1)</li>
198 * <li>fx(1,1)</li>
199 * <li>fy(0,0)</li>
200 * <li>fy(1,0)</li>
201 * <li>fy(0,1)</li>
202 * <li>fy(1,1)</li>
203 * <li>fxy(0,0)</li>
204 * <li>fxy(1,0)</li>
205 * <li>fxy(0,1)</li>
206 * <li>fxy(1,1)</li>
207 * </ul>
208 * @param beta List of function values and function partial derivatives
209 * values
210 * @return the spline coefficients
211 */
212 private double[] computeSplineCoefficients(double[] beta) {
213 final double[] a = new double[16];
214
215 for (int i = 0; i < 16; i++) {
216 double result = 0;
217 final double[] row = AINV[i];
218 for (int j = 0; j < 16; j++) {
219 result += row[j] * beta[j];
220 }
221 a[i] = result;
222 }
223
224 return a;
225 }
226 }
227 /**
228 * 2D-spline function.
229 *
230 * @version $Revision$ $Date$
231 */
232 class BicubicSplineFunction
233 implements BivariateRealFunction {
234 //CHECKSTYLE: stop MultipleVariableDeclarations
235 /** Coefficients */
236 private final double
237 a00, a01, a02, a03,
238 a10, a11, a12, a13,
239 a20, a21, a22, a23,
240 a30, a31, a32, a33;
241 //CHECKSTYLE: resume MultipleVariableDeclarations
242
243 /**
244 * @param a Spline coefficients
245 */
246 public BicubicSplineFunction(double[] a) {
247 this.a00 = a[0];
248 this.a10 = a[1];
249 this.a20 = a[2];
250 this.a30 = a[3];
251 this.a01 = a[4];
252 this.a11 = a[5];
253 this.a21 = a[6];
254 this.a31 = a[7];
255 this.a02 = a[8];
256 this.a12 = a[9];
257 this.a22 = a[10];
258 this.a32 = a[11];
259 this.a03 = a[12];
260 this.a13 = a[13];
261 this.a23 = a[14];
262 this.a33 = a[15];
263 }
264
265 /**
266 * @param x x-coordinate of the interpolation point
267 * @param y y-coordinate of the interpolation point
268 * @return the interpolated value.
269 */
270 public double value(double x, double y) {
271 if (x < 0 || x > 1) {
272 throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
273 x, 0, 1);
274 }
275 if (y < 0 || y > 1) {
276 throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
277 y, 0, 1);
278 }
279
280 final double x2 = x * x;
281 final double x3 = x2 * x;
282 final double y2 = y * y;
283 final double y3 = y2 * y;
284
285 return a00 + a01 * y + a02 * y2 + a03 * y3 +
286 a10 * x + a11 * x * y + a12 * x * y2 + a13 * x * y3 +
287 a20 * x2 + a21 * x2 * y + a22 * x2 * y2 + a23 * x2 * y3 +
288 a30 * x3 + a31 * x3 * y + a32 * x3 * y2 + a33 * x3 * y3;
289 }
290 }