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 org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.Field;
21 import org.hipparchus.analysis.CalculusFieldBivariateFunction;
22 import org.hipparchus.exception.MathIllegalArgumentException;
23 import org.hipparchus.util.MathArrays;
24
25 /**
26 * Interpolate grid data using bi-linear interpolation.
27 * <p>
28 * This interpolator is thread-safe.
29 * </p>
30 * @param <T> Type of the field elements.
31 * @since 4.0
32 */
33 public class FieldBilinearInterpolatingFunction<T extends CalculusFieldElement<T>>
34 implements CalculusFieldBivariateFunction<T> {
35
36 /** Grid along the x axis. */
37 private final FieldGridAxis<T> xGrid;
38
39 /** Grid along the y axis. */
40 private final FieldGridAxis<T> yGrid;
41
42 /** Grid size along the y axis. */
43 private final int ySize;
44
45 /** Values of the interpolation points on all the grid knots (in a flatten array). */
46 private final T[] fVal;
47
48 /** Simple constructor.
49 * @param xVal All the x-coordinates of the interpolation points, sorted
50 * in increasing order.
51 * @param yVal All the y-coordinates of the interpolation points, sorted
52 * in increasing order.
53 * @param fVal The values of the interpolation points on all the grid knots:
54 * {@code fVal[i][j] = f(xVal[i], yVal[j])}.
55 * @exception MathIllegalArgumentException if grid size is smaller than 2
56 * or if the grid is not sorted in strict increasing order
57 */
58 public FieldBilinearInterpolatingFunction(final T[] xVal, final T[] yVal, final T[][] fVal)
59 throws MathIllegalArgumentException {
60 final Field<T> field = fVal[0][0].getField();
61 this.xGrid = new FieldGridAxis<>(xVal, 2);
62 this.yGrid = new FieldGridAxis<>(yVal, 2);
63 this.ySize = yVal.length;
64 this.fVal = MathArrays.buildArray(field, xVal.length * ySize);
65 int k = 0;
66 for (int i = 0; i < xVal.length; ++i) {
67 final T[] fi = fVal[i];
68 for (int j = 0; j < ySize; ++j) {
69 this.fVal[k++] = fi[j];
70 }
71 }
72 }
73
74 /** Get the lowest grid x coordinate.
75 * @return lowest grid x coordinate
76 */
77 public T getXInf() {
78 return xGrid.node(0);
79 }
80
81 /** Get the highest grid x coordinate.
82 * @return highest grid x coordinate
83 */
84 public T getXSup() {
85 return xGrid.node(xGrid.size() - 1);
86 }
87
88 /** Get the lowest grid y coordinate.
89 * @return lowest grid y coordinate
90 */
91 public T getYInf() {
92 return yGrid.node(0);
93 }
94
95 /** Get the highest grid y coordinate.
96 * @return highest grid y coordinate
97 */
98 public T getYSup() {
99 return yGrid.node(yGrid.size() - 1);
100 }
101
102 /** {@inheritDoc} */
103 @Override
104 public T value(T x, T y) {
105
106 // get the interpolation nodes
107 final int i = xGrid.interpolationIndex(x);
108 final int j = yGrid.interpolationIndex(y);
109 final T x0 = xGrid.node(i);
110 final T x1 = xGrid.node(i + 1);
111 final T y0 = yGrid.node(j);
112 final T y1 = yGrid.node(j + 1);
113
114 // get the function values at interpolation nodes
115 final int k0 = i * ySize + j;
116 final int k1 = k0 + ySize;
117 final T z00 = fVal[k0];
118 final T z01 = fVal[k0 + 1];
119 final T z10 = fVal[k1];
120 final T z11 = fVal[k1 + 1];
121
122 // interpolate
123 final T dx0 = x.subtract(x0);
124 final T mdx1 = x.subtract(x1);
125 final T dx10 = x1.subtract(x0);
126 final T dy0 = y.subtract(y0);
127 final T mdy1 = y.subtract(y1);
128 final T dy10 = y1.subtract(y0);
129
130 return dy0.multiply(z11).subtract(mdy1.multiply(z10)).multiply(dx0).
131 subtract(dy0.multiply(z01).subtract(mdy1.multiply(z00)).multiply(mdx1)).
132 divide(dx10.multiply(dy10));
133
134 }
135
136 }