1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 package org.hipparchus.analysis.interpolation;
23
24 import java.util.Arrays;
25
26 import org.hipparchus.CalculusFieldElement;
27 import org.hipparchus.analysis.BivariateFunction;
28 import org.hipparchus.analysis.FieldBivariateFunction;
29 import org.hipparchus.analysis.polynomials.FieldPolynomialSplineFunction;
30 import org.hipparchus.analysis.polynomials.PolynomialSplineFunction;
31 import org.hipparchus.exception.LocalizedCoreFormats;
32 import org.hipparchus.exception.MathIllegalArgumentException;
33 import org.hipparchus.exception.NullArgumentException;
34 import org.hipparchus.util.MathArrays;
35
36
37
38
39
40
41
42
43
44
45
46 public class PiecewiseBicubicSplineInterpolatingFunction
47 implements BivariateFunction, FieldBivariateFunction {
48
49
50 private static final int MIN_NUM_POINTS = 5;
51
52 private final double[] xval;
53
54 private final double[] yval;
55
56 private final double[][] fval;
57
58
59
60
61
62
63
64
65
66
67
68
69
70 public PiecewiseBicubicSplineInterpolatingFunction(double[] x,
71 double[] y,
72 double[][] f)
73 throws MathIllegalArgumentException, NullArgumentException {
74 if (x == null ||
75 y == null ||
76 f == null ||
77 f[0] == null) {
78 throw new NullArgumentException();
79 }
80
81 final int xLen = x.length;
82 final int yLen = y.length;
83
84 if (xLen == 0 ||
85 yLen == 0 ||
86 f.length == 0 ||
87 f[0].length == 0) {
88 throw new MathIllegalArgumentException(LocalizedCoreFormats.NO_DATA);
89 }
90
91 if (xLen < MIN_NUM_POINTS ||
92 yLen < MIN_NUM_POINTS ||
93 f.length < MIN_NUM_POINTS ||
94 f[0].length < MIN_NUM_POINTS) {
95 throw new MathIllegalArgumentException(LocalizedCoreFormats.INSUFFICIENT_DATA);
96 }
97
98 if (xLen != f.length) {
99 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
100 xLen, f.length);
101 }
102
103 if (yLen != f[0].length) {
104 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
105 yLen, f[0].length);
106 }
107
108 MathArrays.checkOrder(x);
109 MathArrays.checkOrder(y);
110
111 xval = x.clone();
112 yval = y.clone();
113 fval = f.clone();
114 }
115
116
117
118
119 @Override
120 public double value(double x,
121 double y)
122 throws MathIllegalArgumentException {
123 final AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator();
124 final int offset = 2;
125 final int count = offset + 3;
126 final int i = searchIndex(x, xval, offset, count);
127 final int j = searchIndex(y, yval, offset, count);
128
129 final double xArray[] = new double[count];
130 final double yArray[] = new double[count];
131 final double zArray[] = new double[count];
132 final double interpArray[] = new double[count];
133
134 for (int index = 0; index < count; index++) {
135 xArray[index] = xval[i + index];
136 yArray[index] = yval[j + index];
137 }
138
139 for (int zIndex = 0; zIndex < count; zIndex++) {
140 for (int index = 0; index < count; index++) {
141 zArray[index] = fval[i + index][j + zIndex];
142 }
143 final PolynomialSplineFunction spline = interpolator.interpolate(xArray, zArray);
144 interpArray[zIndex] = spline.value(x);
145 }
146
147 final PolynomialSplineFunction spline = interpolator.interpolate(yArray, interpArray);
148
149 return spline.value(y);
150
151 }
152
153
154
155
156
157 @Override
158 public <T extends CalculusFieldElement<T>> T value(final T x, final T y)
159 throws MathIllegalArgumentException {
160 final AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator();
161 final int offset = 2;
162 final int count = offset + 3;
163 final int i = searchIndex(x.getReal(), xval, offset, count);
164 final int j = searchIndex(y.getReal(), yval, offset, count);
165
166 final double xArray[] = new double[count];
167 final T yArray[] = MathArrays.buildArray(x.getField(), count);
168 final double zArray[] = new double[count];
169 final T interpArray[] = MathArrays.buildArray(x.getField(), count);
170
171 final T zero = x.getField().getZero();
172 for (int index = 0; index < count; index++) {
173 xArray[index] = xval[i + index];
174 yArray[index] = zero.add(yval[j + index]);
175 }
176
177 for (int zIndex = 0; zIndex < count; zIndex++) {
178 for (int index = 0; index < count; index++) {
179 zArray[index] = fval[i + index][j + zIndex];
180 }
181 final PolynomialSplineFunction spline = interpolator.interpolate(xArray, zArray);
182 interpArray[zIndex] = spline.value(x);
183 }
184
185 final FieldPolynomialSplineFunction<T> spline = interpolator.interpolate(yArray, interpArray);
186
187 return spline.value(y);
188
189 }
190
191
192
193
194
195
196
197
198 public boolean isValidPoint(double x,
199 double y) {
200 if (x < xval[0] ||
201 x > xval[xval.length - 1] ||
202 y < yval[0] ||
203 y > yval[yval.length - 1]) {
204 return false;
205 } else {
206 return true;
207 }
208 }
209
210
211
212
213
214
215
216
217
218
219
220
221 private int searchIndex(double c,
222 double[] val,
223 int offset,
224 int count) {
225 int r = Arrays.binarySearch(val, c);
226
227 if (r == -1 || r == -val.length - 1) {
228 throw new MathIllegalArgumentException(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE,
229 c, val[0], val[val.length - 1]);
230 }
231
232 if (r < 0) {
233
234
235
236 r = -r - offset - 1;
237 } else {
238 r -= offset;
239 }
240
241 if (r < 0) {
242 r = 0;
243 }
244
245 if ((r + count) >= val.length) {
246
247
248 r = val.length - count;
249 }
250
251 return r;
252 }
253 }