| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| package org.apache.commons.math4.legacy.analysis.interpolation; |
| |
| import java.util.Arrays; |
| |
| import org.apache.commons.math4.legacy.analysis.BivariateFunction; |
| import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialSplineFunction; |
| import org.apache.commons.math4.legacy.exception.DimensionMismatchException; |
| import org.apache.commons.math4.legacy.exception.InsufficientDataException; |
| import org.apache.commons.math4.legacy.exception.NoDataException; |
| import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException; |
| import org.apache.commons.math4.legacy.exception.NullArgumentException; |
| import org.apache.commons.math4.legacy.exception.OutOfRangeException; |
| import org.apache.commons.math4.legacy.core.MathArrays; |
| |
| /** |
| * Function that implements the |
| * <a href="http://www.paulinternet.nl/?page=bicubic">bicubic spline</a> |
| * interpolation. |
| * This implementation currently uses {@link AkimaSplineInterpolator} as the |
| * underlying one-dimensional interpolator, which requires 5 sample points; |
| * insufficient data will raise an exception when the |
| * {@link #value(double,double) value} method is called. |
| * |
| * @since 3.4 |
| */ |
| public class PiecewiseBicubicSplineInterpolatingFunction |
| implements BivariateFunction { |
| /** The minimum number of points that are needed to compute the function. */ |
| private static final int MIN_NUM_POINTS = 5; |
| /** Samples x-coordinates. */ |
| private final double[] xval; |
| /** Samples y-coordinates. */ |
| private final double[] yval; |
| /** Set of cubic splines patching the whole data grid. */ |
| private final double[][] fval; |
| |
| /** |
| * @param x Sample values of the x-coordinate, in increasing order. |
| * @param y Sample values of the y-coordinate, in increasing order. |
| * @param f Values of the function on every grid point. the expected number |
| * of elements. |
| * @throws NonMonotonicSequenceException if {@code x} or {@code y} are not |
| * strictly increasing. |
| * @throws NullArgumentException if any of the arguments are null |
| * @throws NoDataException if any of the arrays has zero length. |
| * @throws DimensionMismatchException if the length of x and y don't match the row, column |
| * height of f |
| */ |
| public PiecewiseBicubicSplineInterpolatingFunction(double[] x, |
| double[] y, |
| double[][] f) |
| throws DimensionMismatchException, |
| NullArgumentException, |
| NoDataException, |
| NonMonotonicSequenceException { |
| if (x == null || |
| y == null || |
| f == null || |
| f[0] == null) { |
| throw new NullArgumentException(); |
| } |
| |
| final int xLen = x.length; |
| final int yLen = y.length; |
| |
| if (xLen == 0 || |
| yLen == 0 || |
| f.length == 0 || |
| f[0].length == 0) { |
| throw new NoDataException(); |
| } |
| |
| if (xLen < MIN_NUM_POINTS || |
| yLen < MIN_NUM_POINTS || |
| f.length < MIN_NUM_POINTS || |
| f[0].length < MIN_NUM_POINTS) { |
| throw new InsufficientDataException(); |
| } |
| |
| if (xLen != f.length) { |
| throw new DimensionMismatchException(xLen, f.length); |
| } |
| |
| if (yLen != f[0].length) { |
| throw new DimensionMismatchException(yLen, f[0].length); |
| } |
| |
| MathArrays.checkOrder(x); |
| MathArrays.checkOrder(y); |
| |
| xval = x.clone(); |
| yval = y.clone(); |
| fval = f.clone(); |
| } |
| |
| /** |
| * {@inheritDoc} |
| */ |
| @Override |
| public double value(double x, |
| double y) |
| throws OutOfRangeException { |
| final AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator(); |
| final int offset = 2; |
| final int count = offset + 3; |
| final int i = searchIndex(x, xval, offset, count); |
| final int j = searchIndex(y, yval, offset, count); |
| |
| final double xArray[] = new double[count]; |
| final double yArray[] = new double[count]; |
| final double zArray[] = new double[count]; |
| final double interpArray[] = new double[count]; |
| |
| for (int index = 0; index < count; index++) { |
| xArray[index] = xval[i + index]; |
| yArray[index] = yval[j + index]; |
| } |
| |
| for (int zIndex = 0; zIndex < count; zIndex++) { |
| for (int index = 0; index < count; index++) { |
| zArray[index] = fval[i + index][j + zIndex]; |
| } |
| final PolynomialSplineFunction spline = interpolator.interpolate(xArray, zArray); |
| interpArray[zIndex] = spline.value(x); |
| } |
| |
| final PolynomialSplineFunction spline = interpolator.interpolate(yArray, interpArray); |
| |
| return spline.value(y); |
| } |
| |
| /** |
| * Indicates whether a point is within the interpolation range. |
| * |
| * @param x First coordinate. |
| * @param y Second coordinate. |
| * @return {@code true} if (x, y) is a valid point. |
| * @since 3.3 |
| */ |
| public boolean isValidPoint(double x, |
| double y) { |
| return !(x < xval[0] || |
| x > xval[xval.length - 1] || |
| y < yval[0] || |
| y > yval[yval.length - 1]); |
| } |
| |
| /** |
| * @param c Coordinate. |
| * @param val Coordinate samples. |
| * @param offset how far back from found value to offset for querying |
| * @param count total number of elements forward from beginning that will be |
| * queried |
| * @return the index in {@code val} corresponding to the interval containing |
| * {@code c}. |
| * @throws OutOfRangeException if {@code c} is out of the range defined by |
| * the boundary values of {@code val}. |
| */ |
| private int searchIndex(double c, |
| double[] val, |
| int offset, |
| int count) { |
| int r = Arrays.binarySearch(val, c); |
| |
| if (r == -1 || r == -val.length - 1) { |
| throw new OutOfRangeException(c, val[0], val[val.length - 1]); |
| } |
| |
| if (r < 0) { |
| // "c" in within an interpolation sub-interval, which returns |
| // negative |
| // need to remove the negative sign for consistency |
| r = -r - offset - 1; |
| } else { |
| r -= offset; |
| } |
| |
| if (r < 0) { |
| r = 0; |
| } |
| |
| if ((r + count) >= val.length) { |
| // "c" is the last sample of the range: Return the index |
| // of the sample at the lower end of the last sub-interval. |
| r = val.length - count; |
| } |
| |
| return r; |
| } |
| } |