View Javadoc
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 }