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.analysis.polynomials.PolynomialSplineFunction;
022import org.apache.commons.math3.exception.DimensionMismatchException;
023import org.apache.commons.math3.exception.InsufficientDataException;
024import org.apache.commons.math3.exception.NoDataException;
025import org.apache.commons.math3.exception.NullArgumentException;
026import org.apache.commons.math3.exception.OutOfRangeException;
027import org.apache.commons.math3.exception.NonMonotonicSequenceException;
028import org.apache.commons.math3.util.MathArrays;
029
030/**
031 * Function that implements the <a
032 * href="http://www.paulinternet.nl/?page=bicubic"> bicubic spline
033 * interpolation</a>.
034 *
035 * @since 2.1
036 */
037public class BicubicSplineInterpolatingFunction
038    implements BivariateFunction {
039
040    /** Samples x-coordinates */
041    private final double[] xval;
042
043    /** Samples y-coordinates */
044    private final double[] yval;
045
046    /** Set of cubic splines patching the whole data grid */
047    private final double[][] fval;
048
049    /**
050     * @param x Sample values of the x-coordinate, in increasing order.
051     * @param y Sample values of the y-coordinate, in increasing order.
052     * @param f Values of the function on every grid point. the expected number
053     *        of elements.
054     * @throws NonMonotonicSequenceException if {@code x} or {@code y} are not
055     *         strictly increasing.
056     * @throws NullArgumentException if any of the arguments are null
057     * @throws NoDataException if any of the arrays has zero length.
058     * @throws DimensionMismatchException if the length of x and y don't match the row, column
059     *         height of f
060     */
061
062    public BicubicSplineInterpolatingFunction(double[] x, double[] y,
063                                              double[][] f)
064        throws DimensionMismatchException, NullArgumentException,
065        NoDataException, NonMonotonicSequenceException {
066
067        final int minimumLength = AkimaSplineInterpolator.MINIMUM_NUMBER_POINTS;
068
069        if (x == null || y == null || f == null || f[0] == null) {
070            throw new NullArgumentException();
071        }
072
073        final int xLen = x.length;
074        final int yLen = y.length;
075
076        if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
077            throw new NoDataException();
078        }
079
080        if (xLen < minimumLength || yLen < minimumLength ||
081            f.length < minimumLength || f[0].length < minimumLength) {
082            throw new InsufficientDataException();
083        }
084
085        if (xLen != f.length) {
086            throw new DimensionMismatchException(xLen, f.length);
087        }
088
089        if (yLen != f[0].length) {
090            throw new DimensionMismatchException(yLen, f[0].length);
091        }
092
093        MathArrays.checkOrder(x);
094        MathArrays.checkOrder(y);
095
096        xval = x.clone();
097        yval = y.clone();
098        fval = f.clone();
099    }
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}