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    *      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   */
17  package org.apache.commons.math3.analysis.interpolation;
18  
19  import java.util.Arrays;
20  import org.apache.commons.math3.analysis.BivariateFunction;
21  import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
22  import org.apache.commons.math3.exception.DimensionMismatchException;
23  import org.apache.commons.math3.exception.InsufficientDataException;
24  import org.apache.commons.math3.exception.NoDataException;
25  import org.apache.commons.math3.exception.NullArgumentException;
26  import org.apache.commons.math3.exception.OutOfRangeException;
27  import org.apache.commons.math3.exception.NonMonotonicSequenceException;
28  import org.apache.commons.math3.util.MathArrays;
29  
30  /**
31   * Function that implements the <a
32   * href="http://www.paulinternet.nl/?page=bicubic"> bicubic spline
33   * interpolation</a>.
34   *
35   * @since 2.1
36   */
37  public class BicubicSplineInterpolatingFunction
38      implements BivariateFunction {
39  
40      /** Samples x-coordinates */
41      private final double[] xval;
42  
43      /** Samples y-coordinates */
44      private final double[] yval;
45  
46      /** Set of cubic splines patching the whole data grid */
47      private final double[][] fval;
48  
49      /**
50       * @param x Sample values of the x-coordinate, in increasing order.
51       * @param y Sample values of the y-coordinate, in increasing order.
52       * @param f Values of the function on every grid point. the expected number
53       *        of elements.
54       * @throws NonMonotonicSequenceException if {@code x} or {@code y} are not
55       *         strictly increasing.
56       * @throws NullArgumentException if any of the arguments are null
57       * @throws NoDataException if any of the arrays has zero length.
58       * @throws DimensionMismatchException if the length of x and y don't match the row, column
59       *         height of f
60       */
61  
62      public BicubicSplineInterpolatingFunction(double[] x, double[] y,
63                                                double[][] f)
64          throws DimensionMismatchException, NullArgumentException,
65          NoDataException, NonMonotonicSequenceException {
66  
67          final int minimumLength = AkimaSplineInterpolator.MINIMUM_NUMBER_POINTS;
68  
69          if (x == null || y == null || f == null || f[0] == null) {
70              throw new NullArgumentException();
71          }
72  
73          final int xLen = x.length;
74          final int yLen = y.length;
75  
76          if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
77              throw new NoDataException();
78          }
79  
80          if (xLen < minimumLength || yLen < minimumLength ||
81              f.length < minimumLength || f[0].length < minimumLength) {
82              throw new InsufficientDataException();
83          }
84  
85          if (xLen != f.length) {
86              throw new DimensionMismatchException(xLen, f.length);
87          }
88  
89          if (yLen != f[0].length) {
90              throw new DimensionMismatchException(yLen, f[0].length);
91          }
92  
93          MathArrays.checkOrder(x);
94          MathArrays.checkOrder(y);
95  
96          xval = x.clone();
97          yval = y.clone();
98          fval = f.clone();
99      }
100 
101     /**
102      * {@inheritDoc}
103      */
104     public double value(double x, double y)
105         throws OutOfRangeException {
106         int index;
107         PolynomialSplineFunction spline;
108         AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator();
109         final int offset = 2;
110         final int count = offset + 3;
111         final int i = searchIndex(x, xval, offset, count);
112         final int j = searchIndex(y, yval, offset, count);
113 
114         double xArray[] = new double[count];
115         double yArray[] = new double[count];
116         double zArray[] = new double[count];
117         double interpArray[] = new double[count];
118 
119         for (index = 0; index < count; index++) {
120             xArray[index] = xval[i + index];
121             yArray[index] = yval[j + index];
122         }
123 
124         for (int zIndex = 0; zIndex < count; zIndex++) {
125             for (index = 0; index < count; index++) {
126                 zArray[index] = fval[i + index][j + zIndex];
127             }
128             spline = interpolator.interpolate(xArray, zArray);
129             interpArray[zIndex] = spline.value(x);
130         }
131 
132         spline = interpolator.interpolate(yArray, interpArray);
133 
134         double returnValue = spline.value(y);
135 
136         return returnValue;
137     }
138 
139     /**
140      * Indicates whether a point is within the interpolation range.
141      *
142      * @param x First coordinate.
143      * @param y Second coordinate.
144      * @return {@code true} if (x, y) is a valid point.
145      * @since 3.3
146      */
147     public boolean isValidPoint(double x, double y) {
148         if (x < xval[0] || x > xval[xval.length - 1] || y < yval[0] ||
149             y > yval[yval.length - 1]) {
150             return false;
151         } else {
152             return true;
153         }
154     }
155 
156     /**
157      * @param c Coordinate.
158      * @param val Coordinate samples.
159      * @param offset how far back from found value to offset for querying
160      * @param count total number of elements forward from beginning that will be
161      *        queried
162      * @return the index in {@code val} corresponding to the interval containing
163      *         {@code c}.
164      * @throws OutOfRangeException if {@code c} is out of the range defined by
165      *         the boundary values of {@code val}.
166      */
167     private int searchIndex(double c, double[] val, int offset, int count) {
168         int r = Arrays.binarySearch(val, c);
169 
170         if (r == -1 || r == -val.length - 1) {
171             throw new OutOfRangeException(c, val[0], val[val.length - 1]);
172         }
173 
174         if (r < 0) {
175             // "c" in within an interpolation sub-interval, which returns
176             // negative
177             // need to remove the negative sign for consistency
178             r = -r - offset - 1;
179         } else {
180             r -= offset;
181         }
182 
183         if (r < 0) {
184             r = 0;
185         }
186 
187         if ((r + count) >= val.length) {
188             // "c" is the last sample of the range: Return the index
189             // of the sample at the lower end of the last sub-interval.
190             r = val.length - count;
191         }
192 
193         return r;
194     }
195 }