View Javadoc
1   /*
2    * Licensed to the Hipparchus project 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 Hipparchus project 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  package org.hipparchus.analysis.interpolation;
18  
19  import java.io.Serializable;
20  
21  import org.hipparchus.CalculusFieldElement;
22  import org.hipparchus.analysis.BivariateFunction;
23  import org.hipparchus.analysis.FieldBivariateFunction;
24  import org.hipparchus.exception.MathIllegalArgumentException;
25  
26  /**
27   * Interpolate grid data using bi-linear interpolation.
28   * <p>
29   * This interpolator is thread-safe.
30   * </p>
31   * @since 1.4
32   */
33  public class BilinearInterpolatingFunction implements BivariateFunction, FieldBivariateFunction, Serializable {
34  
35      /** Serializable UID. */
36      private static final long serialVersionUID = 20180926L;
37  
38      /** Grid along the x axis. */
39      private final GridAxis xGrid;
40  
41      /** Grid along the y axis. */
42      private final GridAxis yGrid;
43  
44      /** Grid size along the y axis. */
45      private final int ySize;
46  
47      /** Values of the interpolation points on all the grid knots (in a flatten array). */
48      private final double[] fVal;
49  
50      /** Simple constructor.
51       * @param xVal All the x-coordinates of the interpolation points, sorted
52       * in increasing order.
53       * @param yVal All the y-coordinates of the interpolation points, sorted
54       * in increasing order.
55       * @param fVal The values of the interpolation points on all the grid knots:
56       * {@code fVal[i][j] = f(xVal[i], yVal[j])}.
57       * @exception MathIllegalArgumentException if grid size is smaller than 2
58       * or if the grid is not sorted in strict increasing order
59       */
60      public BilinearInterpolatingFunction(final double[] xVal, final double[] yVal,
61                                           final double[][] fVal)
62          throws MathIllegalArgumentException {
63          this.xGrid = new GridAxis(xVal, 2);
64          this.yGrid = new GridAxis(yVal, 2);
65          this.ySize = yVal.length;
66          this.fVal  = new double[xVal.length * ySize];
67          int k = 0;
68          for (int i = 0; i < xVal.length; ++i) {
69              final double[] fi = fVal[i];
70              for (int j = 0; j < ySize; ++j) {
71                  this.fVal[k++] = fi[j];
72              }
73          }
74      }
75  
76      /** Get the lowest grid x coordinate.
77       * @return lowest grid x coordinate
78       */
79      public double getXInf() {
80          return xGrid.node(0);
81      }
82  
83      /** Get the highest grid x coordinate.
84       * @return highest grid x coordinate
85       */
86      public double getXSup() {
87          return xGrid.node(xGrid.size() - 1);
88      }
89  
90      /** Get the lowest grid y coordinate.
91       * @return lowest grid y coordinate
92       */
93      public double getYInf() {
94          return yGrid.node(0);
95      }
96  
97      /** Get the highest grid y coordinate.
98       * @return highest grid y coordinate
99       */
100     public double getYSup() {
101         return yGrid.node(yGrid.size() - 1);
102     }
103 
104     /** {@inheritDoc} */
105     @Override
106     public double value(final double x, final double y) {
107 
108         // get the interpolation nodes
109         final int    i   = xGrid.interpolationIndex(x);
110         final int    j   = yGrid.interpolationIndex(y);
111         final double x0  = xGrid.node(i);
112         final double x1  = xGrid.node(i + 1);
113         final double y0  = yGrid.node(j);
114         final double y1  = yGrid.node(j + 1);
115 
116         // get the function values at interpolation nodes
117         final int    k0  = i * ySize + j;
118         final int    k1  = k0 + ySize;
119         final double z00 = fVal[k0];
120         final double z01 = fVal[k0 + 1];
121         final double z10 = fVal[k1];
122         final double z11 = fVal[k1 + 1];
123 
124         // interpolate
125         final double dx0  = x  - x0;
126         final double dx1  = x1 - x;
127         final double dx10 = x1 - x0;
128         final double dy0  = y  - y0;
129         final double dy1  = y1 - y;
130         final double dy10 = y1 - y0;
131         return (dx0 * (dy0 * z11 + dy1 * z10) + dx1 * (dy0 * z01 + dy1 * z00)) /
132                (dx10 * dy10);
133 
134     }
135 
136     /** {@inheritDoc}
137      * @since 1.5
138      */
139     @Override
140     public <T extends CalculusFieldElement<T>> T value(T x, T y) {
141 
142         // get the interpolation nodes
143         final int    i   = xGrid.interpolationIndex(x.getReal());
144         final int    j   = yGrid.interpolationIndex(y.getReal());
145         final double x0  = xGrid.node(i);
146         final double x1  = xGrid.node(i + 1);
147         final double y0  = yGrid.node(j);
148         final double y1  = yGrid.node(j + 1);
149 
150         // get the function values at interpolation nodes
151         final int    k0  = i * ySize + j;
152         final int    k1  = k0 + ySize;
153         final double z00 = fVal[k0];
154         final double z01 = fVal[k0 + 1];
155         final double z10 = fVal[k1];
156         final double z11 = fVal[k1 + 1];
157 
158         // interpolate
159         final T      dx0   = x.subtract(x0);
160         final T      mdx1  = x.subtract(x1);
161         final double dx10  = x1 - x0;
162         final T      dy0   = y.subtract(y0);
163         final T      mdy1  = y.subtract(y1);
164         final double dy10  = y1 - y0;
165         return          dy0.multiply(z11).subtract(mdy1.multiply(z10)).multiply(dx0).
166                subtract(dy0.multiply(z01).subtract(mdy1.multiply(z00)).multiply(mdx1)).
167                divide(dx10 * dy10);
168 
169     }
170 
171 }