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.math4.analysis.BivariateFunction;
022import org.apache.commons.math4.analysis.polynomials.PolynomialSplineFunction;
023import org.apache.commons.math4.exception.DimensionMismatchException;
024import org.apache.commons.math4.exception.InsufficientDataException;
025import org.apache.commons.math4.exception.NoDataException;
026import org.apache.commons.math4.exception.NonMonotonicSequenceException;
027import org.apache.commons.math4.exception.NullArgumentException;
028import org.apache.commons.math4.exception.OutOfRangeException;
029import org.apache.commons.math4.util.MathArrays;
030
031/**
032 * Function that implements the
033 * <a href="http://www.paulinternet.nl/?page=bicubic">bicubic spline</a>
034 * interpolation.
035 * This implementation currently uses {@link AkimaSplineInterpolator} as the
036 * underlying one-dimensional interpolator, which requires 5 sample points;
037 * insufficient data will raise an exception when the
038 * {@link #value(double,double) value} method is called.
039 *
040 * @since 3.4
041 */
042public class PiecewiseBicubicSplineInterpolatingFunction
043    implements BivariateFunction {
044    /** The minimum number of points that are needed to compute the function. */
045    private static final int MIN_NUM_POINTS = 5;
046    /** Samples x-coordinates */
047    private final double[] xval;
048    /** Samples y-coordinates */
049    private final double[] yval;
050    /** Set of cubic splines patching the whole data grid */
051    private final double[][] fval;
052
053    /**
054     * @param x Sample values of the x-coordinate, in increasing order.
055     * @param y Sample values of the y-coordinate, in increasing order.
056     * @param f Values of the function on every grid point. the expected number
057     *        of elements.
058     * @throws NonMonotonicSequenceException if {@code x} or {@code y} are not
059     *         strictly increasing.
060     * @throws NullArgumentException if any of the arguments are null
061     * @throws NoDataException if any of the arrays has zero length.
062     * @throws DimensionMismatchException if the length of x and y don't match the row, column
063     *         height of f
064     */
065    public PiecewiseBicubicSplineInterpolatingFunction(double[] x,
066                                                       double[] y,
067                                                       double[][] f)
068        throws DimensionMismatchException,
069               NullArgumentException,
070               NoDataException,
071               NonMonotonicSequenceException {
072        if (x == null ||
073            y == null ||
074            f == null ||
075            f[0] == null) {
076            throw new NullArgumentException();
077        }
078
079        final int xLen = x.length;
080        final int yLen = y.length;
081
082        if (xLen == 0 ||
083            yLen == 0 ||
084            f.length == 0 ||
085            f[0].length == 0) {
086            throw new NoDataException();
087        }
088
089        if (xLen < MIN_NUM_POINTS ||
090            yLen < MIN_NUM_POINTS ||
091            f.length < MIN_NUM_POINTS ||
092            f[0].length < MIN_NUM_POINTS) {
093            throw new InsufficientDataException();
094        }
095
096        if (xLen != f.length) {
097            throw new DimensionMismatchException(xLen, f.length);
098        }
099
100        if (yLen != f[0].length) {
101            throw new DimensionMismatchException(yLen, f[0].length);
102        }
103
104        MathArrays.checkOrder(x);
105        MathArrays.checkOrder(y);
106
107        xval = x.clone();
108        yval = y.clone();
109        fval = f.clone();
110    }
111
112    /**
113     * {@inheritDoc}
114     */
115    @Override
116    public double value(double x,
117                        double y)
118        throws OutOfRangeException {
119        final AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator();
120        final int offset = 2;
121        final int count = offset + 3;
122        final int i = searchIndex(x, xval, offset, count);
123        final int j = searchIndex(y, yval, offset, count);
124
125        final double xArray[] = new double[count];
126        final double yArray[] = new double[count];
127        final double zArray[] = new double[count];
128        final double interpArray[] = new double[count];
129
130        for (int index = 0; index < count; index++) {
131            xArray[index] = xval[i + index];
132            yArray[index] = yval[j + index];
133        }
134
135        for (int zIndex = 0; zIndex < count; zIndex++) {
136            for (int index = 0; index < count; index++) {
137                zArray[index] = fval[i + index][j + zIndex];
138            }
139            final PolynomialSplineFunction spline = interpolator.interpolate(xArray, zArray);
140            interpArray[zIndex] = spline.value(x);
141        }
142
143        final PolynomialSplineFunction spline = interpolator.interpolate(yArray, interpArray);
144
145        double returnValue = spline.value(y);
146
147        return returnValue;
148    }
149
150    /**
151     * Indicates whether a point is within the interpolation range.
152     *
153     * @param x First coordinate.
154     * @param y Second coordinate.
155     * @return {@code true} if (x, y) is a valid point.
156     * @since 3.3
157     */
158    public boolean isValidPoint(double x,
159                                double y) {
160        if (x < xval[0] ||
161            x > xval[xval.length - 1] ||
162            y < yval[0] ||
163            y > yval[yval.length - 1]) {
164            return false;
165        } else {
166            return true;
167        }
168    }
169
170    /**
171     * @param c Coordinate.
172     * @param val Coordinate samples.
173     * @param offset how far back from found value to offset for querying
174     * @param count total number of elements forward from beginning that will be
175     *        queried
176     * @return the index in {@code val} corresponding to the interval containing
177     *         {@code c}.
178     * @throws OutOfRangeException if {@code c} is out of the range defined by
179     *         the boundary values of {@code val}.
180     */
181    private int searchIndex(double c,
182                            double[] val,
183                            int offset,
184                            int count) {
185        int r = Arrays.binarySearch(val, c);
186
187        if (r == -1 || r == -val.length - 1) {
188            throw new OutOfRangeException(c, val[0], val[val.length - 1]);
189        }
190
191        if (r < 0) {
192            // "c" in within an interpolation sub-interval, which returns
193            // negative
194            // need to remove the negative sign for consistency
195            r = -r - offset - 1;
196        } else {
197            r -= offset;
198        }
199
200        if (r < 0) {
201            r = 0;
202        }
203
204        if ((r + count) >= val.length) {
205            // "c" is the last sample of the range: Return the index
206            // of the sample at the lower end of the last sub-interval.
207            r = val.length - count;
208        }
209
210        return r;
211    }
212}