View Javadoc
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    *      https://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   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.analysis.interpolation;
23  
24  import java.util.Arrays;
25  
26  import org.hipparchus.analysis.BivariateFunction;
27  import org.hipparchus.exception.LocalizedCoreFormats;
28  import org.hipparchus.exception.MathIllegalArgumentException;
29  import org.hipparchus.util.MathArrays;
30  import org.hipparchus.util.MathUtils;
31  
32  /**
33   * Function that implements the
34   * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
35   * bicubic spline interpolation</a>.
36   *
37   */
38  public class BicubicInterpolatingFunction
39      implements BivariateFunction {
40      /** Number of coefficients. */
41      private static final int NUM_COEFF = 16;
42      /**
43       * Matrix to compute the spline coefficients from the function values
44       * and function derivatives values
45       */
46      private static final double[][] AINV = {
47          { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
48          { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
49          { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
50          { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
51          { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
52          { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
53          { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
54          { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
55          { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
56          { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
57          { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
58          { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
59          { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
60          { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
61          { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
62          { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
63      };
64  
65      /** Samples x-coordinates */
66      private final double[] xval;
67      /** Samples y-coordinates */
68      private final double[] yval;
69      /** Set of cubic splines patching the whole data grid */
70      private final BicubicFunction[][] splines;
71  
72      /** Simple constructor.
73       * @param x Sample values of the x-coordinate, in increasing order.
74       * @param y Sample values of the y-coordinate, in increasing order.
75       * @param f Values of the function on every grid point.
76       * @param dFdX Values of the partial derivative of function with respect
77       * to x on every grid point.
78       * @param dFdY Values of the partial derivative of function with respect
79       * to y on every grid point.
80       * @param d2FdXdY Values of the cross partial derivative of function on
81       * every grid point.
82       * @throws MathIllegalArgumentException if the various arrays do not contain
83       * the expected number of elements.
84       * @throws MathIllegalArgumentException if {@code x} or {@code y} are
85       * not strictly increasing.
86       * @throws MathIllegalArgumentException if any of the arrays has zero length.
87       */
88      public BicubicInterpolatingFunction(double[] x,
89                                          double[] y,
90                                          double[][] f,
91                                          double[][] dFdX,
92                                          double[][] dFdY,
93                                          double[][] d2FdXdY)
94          throws MathIllegalArgumentException {
95          final int xLen = x.length;
96          final int yLen = y.length;
97  
98          if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
99              throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
100         }
101         MathUtils.checkDimension(xLen, f.length);
102         MathUtils.checkDimension(xLen, dFdX.length);
103         MathUtils.checkDimension(xLen, dFdY.length);
104         MathUtils.checkDimension(xLen, d2FdXdY.length);
105         MathArrays.checkOrder(x);
106         MathArrays.checkOrder(y);
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 BicubicFunction[lastI][lastJ];
114 
115         for (int i = 0; i < lastI; i++) {
116             MathUtils.checkDimension(f[i].length, yLen);
117             MathUtils.checkDimension(dFdX[i].length, yLen);
118             MathUtils.checkDimension(dFdY[i].length, yLen);
119             MathUtils.checkDimension(d2FdXdY[i].length, yLen);
120 
121             final int ip1 = i + 1;
122             final double xR = xval[ip1] - xval[i];
123             for (int j = 0; j < lastJ; j++) {
124                 final int jp1 = j + 1;
125                 final double yR = yval[jp1] - yval[j];
126                 final double xRyR = xR * yR;
127                 final double[] beta = {
128                     f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
129                     dFdX[i][j] * xR, dFdX[ip1][j] * xR, dFdX[i][jp1] * xR, dFdX[ip1][jp1] * xR,
130                     dFdY[i][j] * yR, dFdY[ip1][j] * yR, dFdY[i][jp1] * yR, dFdY[ip1][jp1] * yR,
131                     d2FdXdY[i][j] * xRyR, d2FdXdY[ip1][j] * xRyR, d2FdXdY[i][jp1] * xRyR, d2FdXdY[ip1][jp1] * xRyR
132                 };
133 
134                 splines[i][j] = new BicubicFunction(computeSplineCoefficients(beta));
135             }
136         }
137     }
138 
139     /**
140      * {@inheritDoc}
141      */
142     @Override
143     public double value(double x, double y)
144         throws MathIllegalArgumentException {
145         final int i = searchIndex(x, xval);
146         final int j = searchIndex(y, yval);
147 
148         final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
149         final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
150 
151         return splines[i][j].value(xN, yN);
152     }
153 
154     /**
155      * Indicates whether a point is within the interpolation range.
156      *
157      * @param x First coordinate.
158      * @param y Second coordinate.
159      * @return {@code true} if (x, y) is a valid point.
160      */
161     public boolean isValidPoint(double x, double y) {
162         if (x < xval[0] ||
163             x > xval[xval.length - 1] ||
164             y < yval[0] ||
165             y > yval[yval.length - 1]) {
166             return false;
167         } else {
168             return true;
169         }
170     }
171 
172     /**
173      * @param c Coordinate.
174      * @param val Coordinate samples.
175      * @return the index in {@code val} corresponding to the interval
176      * containing {@code c}.
177      * @throws MathIllegalArgumentException if {@code c} is out of the
178      * range defined by the boundary values of {@code val}.
179      */
180     private int searchIndex(double c, double[] val) {
181         final int r = Arrays.binarySearch(val, c);
182 
183         if (r == -1 ||
184             r == -val.length - 1) {
185             throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
186                                                    c, val[0], val[val.length - 1]);
187         }
188 
189         if (r < 0) {
190             // "c" in within an interpolation sub-interval: Return the
191             // index of the sample at the lower end of the sub-interval.
192             return -r - 2;
193         }
194         final int last = val.length - 1;
195         if (r == last) {
196             // "c" is the last sample of the range: Return the index
197             // of the sample at the lower end of the last sub-interval.
198             return last - 1;
199         }
200 
201         // "c" is another sample point.
202         return r;
203     }
204 
205     /**
206      * Compute the spline coefficients from the list of function values and
207      * function partial derivatives values at the four corners of a grid
208      * element. They must be specified in the following order:
209      * <ul>
210      *  <li>f(0,0)</li>
211      *  <li>f(1,0)</li>
212      *  <li>f(0,1)</li>
213      *  <li>f(1,1)</li>
214      *  <li>f<sub>x</sub>(0,0)</li>
215      *  <li>f<sub>x</sub>(1,0)</li>
216      *  <li>f<sub>x</sub>(0,1)</li>
217      *  <li>f<sub>x</sub>(1,1)</li>
218      *  <li>f<sub>y</sub>(0,0)</li>
219      *  <li>f<sub>y</sub>(1,0)</li>
220      *  <li>f<sub>y</sub>(0,1)</li>
221      *  <li>f<sub>y</sub>(1,1)</li>
222      *  <li>f<sub>xy</sub>(0,0)</li>
223      *  <li>f<sub>xy</sub>(1,0)</li>
224      *  <li>f<sub>xy</sub>(0,1)</li>
225      *  <li>f<sub>xy</sub>(1,1)</li>
226      * </ul>
227      * where the subscripts indicate the partial derivative with respect to
228      * the corresponding variable(s).
229      *
230      * @param beta List of function values and function partial derivatives
231      * values.
232      * @return the spline coefficients.
233      */
234     private double[] computeSplineCoefficients(double[] beta) {
235         final double[] a = new double[NUM_COEFF];
236 
237         for (int i = 0; i < NUM_COEFF; i++) {
238             double result = 0;
239             final double[] row = AINV[i];
240             for (int j = 0; j < NUM_COEFF; j++) {
241                 result += row[j] * beta[j];
242             }
243             a[i] = result;
244         }
245 
246         return a;
247     }
248 }
249 
250 /**
251  * Bicubic function.
252  */
253 class BicubicFunction implements BivariateFunction {
254     /** Number of points. */
255     private static final short N = 4;
256     /** Coefficients */
257     private final double[][] a;
258 
259     /**
260      * Simple constructor.
261      *
262      * @param coeff Spline coefficients.
263      */
264     BicubicFunction(double[] coeff) {
265         a = new double[N][N];
266         for (int j = 0; j < N; j++) {
267             final double[] aJ = a[j];
268             for (int i = 0; i < N; i++) {
269                 aJ[i] = coeff[i * N + j];
270             }
271         }
272     }
273 
274     /**
275      * {@inheritDoc}
276      */
277     @Override
278     public double value(double x, double y) {
279         MathUtils.checkRangeInclusive(x, 0, 1);
280         MathUtils.checkRangeInclusive(y, 0, 1);
281 
282         final double x2 = x * x;
283         final double x3 = x2 * x;
284         final double[] pX = {1, x, x2, x3};
285 
286         final double y2 = y * y;
287         final double y3 = y2 * y;
288         final double[] pY = {1, y, y2, y3};
289 
290         return apply(pX, pY, a);
291     }
292 
293     /**
294      * Compute the value of the bicubic polynomial.
295      *
296      * @param pX Powers of the x-coordinate.
297      * @param pY Powers of the y-coordinate.
298      * @param coeff Spline coefficients.
299      * @return the interpolated value.
300      */
301     private double apply(double[] pX, double[] pY, double[][] coeff) {
302         double result = 0;
303         for (int i = 0; i < N; i++) {
304             final double r = MathArrays.linearCombination(coeff[i], pY);
305             result += r * pX[i];
306         }
307 
308         return result;
309     }
310 }