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 */
017package org.apache.commons.math4.analysis.interpolation;
018
019import java.util.Arrays;
020
021import org.apache.commons.numbers.arrays.LinearCombination;
022import org.apache.commons.math4.analysis.BivariateFunction;
023import org.apache.commons.math4.exception.DimensionMismatchException;
024import org.apache.commons.math4.exception.NoDataException;
025import org.apache.commons.math4.exception.NonMonotonicSequenceException;
026import org.apache.commons.math4.exception.OutOfRangeException;
027import org.apache.commons.math4.util.MathArrays;
028
029/**
030 * Function that implements the
031 * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
032 * bicubic spline interpolation</a>.
033 *
034 * @since 3.4
035 */
036public class BicubicInterpolatingFunction
037    implements BivariateFunction {
038    /** Number of coefficients. */
039    private static final int NUM_COEFF = 16;
040    /**
041     * Matrix to compute the spline coefficients from the function values
042     * and function derivatives values
043     */
044    private static final double[][] AINV = {
045        { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
046        { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
047        { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
048        { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
049        { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
050        { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
051        { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
052        { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
053        { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
054        { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
055        { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
056        { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
057        { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
058        { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
059        { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
060        { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
061    };
062
063    /** Samples x-coordinates */
064    private final double[] xval;
065    /** Samples y-coordinates */
066    private final double[] yval;
067    /** Set of cubic splines patching the whole data grid */
068    private final BicubicFunction[][] splines;
069
070    /**
071     * @param x Sample values of the x-coordinate, in increasing order.
072     * @param y Sample values of the y-coordinate, in increasing order.
073     * @param f Values of the function on every grid point.
074     * @param dFdX Values of the partial derivative of function with respect
075     * to x on every grid point.
076     * @param dFdY Values of the partial derivative of function with respect
077     * to y on every grid point.
078     * @param d2FdXdY Values of the cross partial derivative of function on
079     * every grid point.
080     * @throws DimensionMismatchException if the various arrays do not contain
081     * the expected number of elements.
082     * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
083     * not strictly increasing.
084     * @throws NoDataException if any of the arrays has zero length.
085     */
086    public BicubicInterpolatingFunction(double[] x,
087                                        double[] y,
088                                        double[][] f,
089                                        double[][] dFdX,
090                                        double[][] dFdY,
091                                        double[][] d2FdXdY)
092        throws DimensionMismatchException,
093               NoDataException,
094               NonMonotonicSequenceException {
095        final int xLen = x.length;
096        final int yLen = y.length;
097
098        if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
099            throw new NoDataException();
100        }
101        if (xLen != f.length) {
102            throw new DimensionMismatchException(xLen, f.length);
103        }
104        if (xLen != dFdX.length) {
105            throw new DimensionMismatchException(xLen, dFdX.length);
106        }
107        if (xLen != dFdY.length) {
108            throw new DimensionMismatchException(xLen, dFdY.length);
109        }
110        if (xLen != d2FdXdY.length) {
111            throw new DimensionMismatchException(xLen, d2FdXdY.length);
112        }
113
114        MathArrays.checkOrder(x);
115        MathArrays.checkOrder(y);
116
117        xval = x.clone();
118        yval = y.clone();
119
120        final int lastI = xLen - 1;
121        final int lastJ = yLen - 1;
122        splines = new BicubicFunction[lastI][lastJ];
123
124        for (int i = 0; i < lastI; i++) {
125            if (f[i].length != yLen) {
126                throw new DimensionMismatchException(f[i].length, yLen);
127            }
128            if (dFdX[i].length != yLen) {
129                throw new DimensionMismatchException(dFdX[i].length, yLen);
130            }
131            if (dFdY[i].length != yLen) {
132                throw new DimensionMismatchException(dFdY[i].length, yLen);
133            }
134            if (d2FdXdY[i].length != yLen) {
135                throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
136            }
137            final int ip1 = i + 1;
138            final double xR = xval[ip1] - xval[i];
139            for (int j = 0; j < lastJ; j++) {
140                final int jp1 = j + 1;
141                final double yR = yval[jp1] - yval[j];
142                final double xRyR = xR * yR;
143                final double[] beta = new double[] {
144                    f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
145                    dFdX[i][j] * xR, dFdX[ip1][j] * xR, dFdX[i][jp1] * xR, dFdX[ip1][jp1] * xR,
146                    dFdY[i][j] * yR, dFdY[ip1][j] * yR, dFdY[i][jp1] * yR, dFdY[ip1][jp1] * yR,
147                    d2FdXdY[i][j] * xRyR, d2FdXdY[ip1][j] * xRyR, d2FdXdY[i][jp1] * xRyR, d2FdXdY[ip1][jp1] * xRyR
148                };
149
150                splines[i][j] = new BicubicFunction(computeSplineCoefficients(beta));
151            }
152        }
153    }
154
155    /**
156     * {@inheritDoc}
157     */
158    @Override
159    public double value(double x, double y)
160        throws OutOfRangeException {
161        final int i = searchIndex(x, xval);
162        final int j = searchIndex(y, yval);
163
164        final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
165        final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
166
167        return splines[i][j].value(xN, yN);
168    }
169
170    /**
171     * Indicates whether a point is within the interpolation range.
172     *
173     * @param x First coordinate.
174     * @param y Second coordinate.
175     * @return {@code true} if (x, y) is a valid point.
176     */
177    public boolean isValidPoint(double x, double y) {
178        if (x < xval[0] ||
179            x > xval[xval.length - 1] ||
180            y < yval[0] ||
181            y > yval[yval.length - 1]) {
182            return false;
183        } else {
184            return true;
185        }
186    }
187
188    /**
189     * @param c Coordinate.
190     * @param val Coordinate samples.
191     * @return the index in {@code val} corresponding to the interval
192     * containing {@code c}.
193     * @throws OutOfRangeException if {@code c} is out of the
194     * range defined by the boundary values of {@code val}.
195     */
196    private int searchIndex(double c, double[] val) {
197        final int r = Arrays.binarySearch(val, c);
198
199        if (r == -1 ||
200            r == -val.length - 1) {
201            throw new OutOfRangeException(c, val[0], val[val.length - 1]);
202        }
203
204        if (r < 0) {
205            // "c" in within an interpolation sub-interval: Return the
206            // index of the sample at the lower end of the sub-interval.
207            return -r - 2;
208        }
209        final int last = val.length - 1;
210        if (r == last) {
211            // "c" is the last sample of the range: Return the index
212            // of the sample at the lower end of the last sub-interval.
213            return last - 1;
214        }
215
216        // "c" is another sample point.
217        return r;
218    }
219
220    /**
221     * Compute the spline coefficients from the list of function values and
222     * function partial derivatives values at the four corners of a grid
223     * element. They must be specified in the following order:
224     * <ul>
225     *  <li>f(0,0)</li>
226     *  <li>f(1,0)</li>
227     *  <li>f(0,1)</li>
228     *  <li>f(1,1)</li>
229     *  <li>f<sub>x</sub>(0,0)</li>
230     *  <li>f<sub>x</sub>(1,0)</li>
231     *  <li>f<sub>x</sub>(0,1)</li>
232     *  <li>f<sub>x</sub>(1,1)</li>
233     *  <li>f<sub>y</sub>(0,0)</li>
234     *  <li>f<sub>y</sub>(1,0)</li>
235     *  <li>f<sub>y</sub>(0,1)</li>
236     *  <li>f<sub>y</sub>(1,1)</li>
237     *  <li>f<sub>xy</sub>(0,0)</li>
238     *  <li>f<sub>xy</sub>(1,0)</li>
239     *  <li>f<sub>xy</sub>(0,1)</li>
240     *  <li>f<sub>xy</sub>(1,1)</li>
241     * </ul>
242     * where the subscripts indicate the partial derivative with respect to
243     * the corresponding variable(s).
244     *
245     * @param beta List of function values and function partial derivatives
246     * values.
247     * @return the spline coefficients.
248     */
249    private double[] computeSplineCoefficients(double[] beta) {
250        final double[] a = new double[NUM_COEFF];
251
252        for (int i = 0; i < NUM_COEFF; i++) {
253            double result = 0;
254            final double[] row = AINV[i];
255            for (int j = 0; j < NUM_COEFF; j++) {
256                result += row[j] * beta[j];
257            }
258            a[i] = result;
259        }
260
261        return a;
262    }
263}
264
265/**
266 * Bicubic function.
267 */
268class BicubicFunction implements BivariateFunction {
269    /** Number of points. */
270    private static final short N = 4;
271    /** Coefficients */
272    private final double[][] a;
273
274    /**
275     * Simple constructor.
276     *
277     * @param coeff Spline coefficients.
278     */
279    BicubicFunction(double[] coeff) {
280        a = new double[N][N];
281        for (int j = 0; j < N; j++) {
282            final double[] aJ = a[j];
283            for (int i = 0; i < N; i++) {
284                aJ[i] = coeff[i * N + j];
285            }
286        }
287    }
288
289    /**
290     * {@inheritDoc}
291     */
292    @Override
293    public double value(double x, double y) {
294        if (x < 0 || x > 1) {
295            throw new OutOfRangeException(x, 0, 1);
296        }
297        if (y < 0 || y > 1) {
298            throw new OutOfRangeException(y, 0, 1);
299        }
300
301        final double x2 = x * x;
302        final double x3 = x2 * x;
303        final double[] pX = {1, x, x2, x3};
304
305        final double y2 = y * y;
306        final double y3 = y2 * y;
307        final double[] pY = {1, y, y2, y3};
308
309        return apply(pX, pY, a);
310    }
311
312    /**
313     * Compute the value of the bicubic polynomial.
314     *
315     * @param pX Powers of the x-coordinate.
316     * @param pY Powers of the y-coordinate.
317     * @param coeff Spline coefficients.
318     * @return the interpolated value.
319     */
320    private double apply(double[] pX, double[] pY, double[][] coeff) {
321        double result = 0;
322        for (int i = 0; i < N; i++) {
323            final double r = LinearCombination.value(coeff[i], pY);
324            result += r * pX[i];
325        }
326
327        return result;
328    }
329}