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