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.legacy.analysis.interpolation;
018
019import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
020import org.apache.commons.math4.legacy.exception.NoDataException;
021import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
022import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
023import org.apache.commons.math4.legacy.core.MathArrays;
024
025/**
026 * Generates a {@link BicubicInterpolatingFunction bicubic interpolating
027 * function}.
028 * <p>
029 *  Caveat: Because the interpolation scheme requires that derivatives be
030 *  specified at the sample points, those are approximated with finite
031 *  differences (using the 2-points symmetric formulae).
032 *  Since their values are undefined at the borders of the provided
033 *  interpolation ranges, the interpolated values will be wrong at the
034 *  edges of the patch.
035 *  The {@code interpolate} method will return a function that overrides
036 *  {@link BicubicInterpolatingFunction#isValidPoint(double,double)} to
037 *  indicate points where the interpolation will be inaccurate.
038 * </p>
039 *
040 * @since 3.4
041 */
042public class BicubicInterpolator
043    implements BivariateGridInterpolator {
044    /**
045    * Whether to initialize internal data used to compute the analytical
046    * derivatives of the splines.
047    */
048    private final boolean initializeDerivatives;
049
050    /**
051     * Default constructor.
052     * The argument {@link #BicubicInterpolator(boolean) initializeDerivatives}
053     * is set to {@code false}.
054     */
055    public BicubicInterpolator() {
056        this(false);
057    }
058
059    /**
060     * Creates an interpolator.
061     *
062     * @param initializeDerivatives Whether to initialize the internal data
063     * needed for calling any of the methods that compute the partial derivatives
064     * of the {@link BicubicInterpolatingFunction function} returned from
065     * the call to {@link #interpolate(double[],double[],double[][]) interpolate}.
066     */
067    public BicubicInterpolator(boolean initializeDerivatives) {
068        this.initializeDerivatives = initializeDerivatives;
069    }
070    /**
071     * {@inheritDoc}
072     */
073    @Override
074    public BicubicInterpolatingFunction interpolate(final double[] xval,
075                                                    final double[] yval,
076                                                    final double[][] fval)
077        throws NoDataException, DimensionMismatchException,
078               NonMonotonicSequenceException, NumberIsTooSmallException {
079        if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
080            throw new NoDataException();
081        }
082        if (xval.length != fval.length) {
083            throw new DimensionMismatchException(xval.length, fval.length);
084        }
085
086        MathArrays.checkOrder(xval);
087        MathArrays.checkOrder(yval);
088
089        final int xLen = xval.length;
090        final int yLen = yval.length;
091
092        // Approximation to the partial derivatives using finite differences.
093        final double[][] dFdX = new double[xLen][yLen];
094        final double[][] dFdY = new double[xLen][yLen];
095        final double[][] d2FdXdY = new double[xLen][yLen];
096        for (int i = 1; i < xLen - 1; i++) {
097            final int nI = i + 1;
098            final int pI = i - 1;
099
100            final double nX = xval[nI];
101            final double pX = xval[pI];
102
103            final double deltaX = nX - pX;
104
105            for (int j = 1; j < yLen - 1; j++) {
106                final int nJ = j + 1;
107                final int pJ = j - 1;
108
109                final double nY = yval[nJ];
110                final double pY = yval[pJ];
111
112                final double deltaY = nY - pY;
113
114                dFdX[i][j] = (fval[nI][j] - fval[pI][j]) / deltaX;
115                dFdY[i][j] = (fval[i][nJ] - fval[i][pJ]) / deltaY;
116
117                final double deltaXY = deltaX * deltaY;
118
119                d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] - fval[pI][nJ] + fval[pI][pJ]) / deltaXY;
120            }
121        }
122
123        // Create the interpolating function.
124        return new BicubicInterpolatingFunction(xval, yval, fval,
125                                                dFdX, dFdY, d2FdXdY,
126                                                initializeDerivatives) {
127            /** {@inheritDoc} */
128            @Override
129            public boolean isValidPoint(double x, double y) {
130                return !(x < xval[1] ||
131                    x > xval[xval.length - 2] ||
132                    y < yval[1] ||
133                    y > yval[yval.length - 2]);
134            }
135        };
136    }
137}