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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.analysis.interpolation;
23
24 import org.hipparchus.exception.LocalizedCoreFormats;
25 import org.hipparchus.exception.MathIllegalArgumentException;
26 import org.hipparchus.util.MathArrays;
27 import org.hipparchus.util.MathUtils;
28
29 /**
30 * Generates a {@link BicubicInterpolatingFunction bicubic interpolating
31 * function}.
32 * <p>
33 * Caveat: Because the interpolation scheme requires that derivatives be
34 * specified at the sample points, those are approximated with finite
35 * differences (using the 2-points symmetric formulae).
36 * Since their values are undefined at the borders of the provided
37 * interpolation ranges, the interpolated values will be wrong at the
38 * edges of the patch.
39 * The {@code interpolate} method will return a function that overrides
40 * {@link BicubicInterpolatingFunction#isValidPoint(double,double)} to
41 * indicate points where the interpolation will be inaccurate.
42 * </p>
43 *
44 */
45 public class BicubicInterpolator
46 implements BivariateGridInterpolator {
47
48 /** Empty constructor.
49 * <p>
50 * This constructor is not strictly necessary, but it prevents spurious
51 * javadoc warnings with JDK 18 and later.
52 * </p>
53 * @since 3.0
54 */
55 public BicubicInterpolator() { // NOPMD - unnecessary constructor added intentionally to make javadoc happy
56 // nothing to do
57 }
58
59 /**
60 * {@inheritDoc}
61 */
62 @Override
63 public BicubicInterpolatingFunction interpolate(final double[] xval,
64 final double[] yval,
65 final double[][] fval)
66 throws MathIllegalArgumentException {
67 if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
68 throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
69 }
70 MathUtils.checkDimension(xval.length, fval.length);
71 MathArrays.checkOrder(xval);
72 MathArrays.checkOrder(yval);
73
74 final int xLen = xval.length;
75 final int yLen = yval.length;
76
77 // Approximation to the partial derivatives using finite differences.
78 final double[][] dFdX = new double[xLen][yLen];
79 final double[][] dFdY = new double[xLen][yLen];
80 final double[][] d2FdXdY = new double[xLen][yLen];
81 for (int i = 1; i < xLen - 1; i++) {
82 final int nI = i + 1;
83 final int pI = i - 1;
84
85 final double nX = xval[nI];
86 final double pX = xval[pI];
87
88 final double deltaX = nX - pX;
89
90 for (int j = 1; j < yLen - 1; j++) {
91 final int nJ = j + 1;
92 final int pJ = j - 1;
93
94 final double nY = yval[nJ];
95 final double pY = yval[pJ];
96
97 final double deltaY = nY - pY;
98
99 dFdX[i][j] = (fval[nI][j] - fval[pI][j]) / deltaX;
100 dFdY[i][j] = (fval[i][nJ] - fval[i][pJ]) / deltaY;
101
102 final double deltaXY = deltaX * deltaY;
103
104 d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] - fval[pI][nJ] + fval[pI][pJ]) / deltaXY;
105 }
106 }
107
108 // Create the interpolating function.
109 return new BicubicInterpolatingFunction(xval, yval, fval,
110 dFdX, dFdY, d2FdXdY) {
111 /** {@inheritDoc} */
112 @Override
113 public boolean isValidPoint(double x, double y) {
114 if (x < xval[1] ||
115 x > xval[xval.length - 2] ||
116 y < yval[1] ||
117 y > yval[yval.length - 2]) {
118 return false;
119 } else {
120 return true;
121 }
122 }
123 };
124 }
125 }