1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math4.legacy.analysis.interpolation;
18
19 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
20 import org.apache.commons.math4.legacy.exception.NoDataException;
21 import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
22 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
23 import org.apache.commons.math4.legacy.core.MathArrays;
24
25 /**
26 * Generates a {@link BicubicInterpolatingFunction bicubic interpolating
27 * function}.
28 * <p>
29 * Caveat: Because the interpolation scheme requires that derivatives be
30 * specified at the sample points, those are approximated with finite
31 * differences (using the 2-points symmetric formulae).
32 * Since their values are undefined at the borders of the provided
33 * interpolation ranges, the interpolated values will be wrong at the
34 * edges of the patch.
35 * The {@code interpolate} method will return a function that overrides
36 * {@link BicubicInterpolatingFunction#isValidPoint(double,double)} to
37 * indicate points where the interpolation will be inaccurate.
38 * </p>
39 *
40 * @since 3.4
41 */
42 public class BicubicInterpolator
43 implements BivariateGridInterpolator {
44 /**
45 * Whether to initialize internal data used to compute the analytical
46 * derivatives of the splines.
47 */
48 private final boolean initializeDerivatives;
49
50 /**
51 * Default constructor.
52 * The argument {@link #BicubicInterpolator(boolean) initializeDerivatives}
53 * is set to {@code false}.
54 */
55 public BicubicInterpolator() {
56 this(false);
57 }
58
59 /**
60 * Creates an interpolator.
61 *
62 * @param initializeDerivatives Whether to initialize the internal data
63 * needed for calling any of the methods that compute the partial derivatives
64 * of the {@link BicubicInterpolatingFunction function} returned from
65 * the call to {@link #interpolate(double[],double[],double[][]) interpolate}.
66 */
67 public BicubicInterpolator(boolean initializeDerivatives) {
68 this.initializeDerivatives = initializeDerivatives;
69 }
70 /**
71 * {@inheritDoc}
72 */
73 @Override
74 public BicubicInterpolatingFunction interpolate(final double[] xval,
75 final double[] yval,
76 final double[][] fval)
77 throws NoDataException, DimensionMismatchException,
78 NonMonotonicSequenceException, NumberIsTooSmallException {
79 if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
80 throw new NoDataException();
81 }
82 if (xval.length != fval.length) {
83 throw new DimensionMismatchException(xval.length, fval.length);
84 }
85
86 MathArrays.checkOrder(xval);
87 MathArrays.checkOrder(yval);
88
89 final int xLen = xval.length;
90 final int yLen = yval.length;
91
92 // Approximation to the partial derivatives using finite differences.
93 final double[][] dFdX = new double[xLen][yLen];
94 final double[][] dFdY = new double[xLen][yLen];
95 final double[][] d2FdXdY = new double[xLen][yLen];
96 for (int i = 1; i < xLen - 1; i++) {
97 final int nI = i + 1;
98 final int pI = i - 1;
99
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 }