commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunction.java

/*
 * 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;
    }
}