1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math4.legacy.analysis.interpolation;
18
19 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
20 import org.apache.commons.math4.legacy.exception.NoDataException;
21 import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
22 import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
23 import org.apache.commons.math4.legacy.core.MathArrays;
24
25
26
27
28
29
30 public class TricubicInterpolator
31 implements TrivariateGridInterpolator {
32
33
34
35 @Override
36 public TricubicInterpolatingFunction interpolate(final double[] xval,
37 final double[] yval,
38 final double[] zval,
39 final double[][][] fval)
40 throws NoDataException, NumberIsTooSmallException,
41 DimensionMismatchException, NonMonotonicSequenceException {
42 if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) {
43 throw new NoDataException();
44 }
45 if (xval.length != fval.length) {
46 throw new DimensionMismatchException(xval.length, fval.length);
47 }
48
49 MathArrays.checkOrder(xval);
50 MathArrays.checkOrder(yval);
51 MathArrays.checkOrder(zval);
52
53 final int xLen = xval.length;
54 final int yLen = yval.length;
55 final int zLen = zval.length;
56
57
58 final double[][][] dFdX = new double[xLen][yLen][zLen];
59 final double[][][] dFdY = new double[xLen][yLen][zLen];
60 final double[][][] dFdZ = new double[xLen][yLen][zLen];
61 final double[][][] d2FdXdY = new double[xLen][yLen][zLen];
62 final double[][][] d2FdXdZ = new double[xLen][yLen][zLen];
63 final double[][][] d2FdYdZ = new double[xLen][yLen][zLen];
64 final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen];
65
66 for (int i = 1; i < xLen - 1; i++) {
67 if (yval.length != fval[i].length) {
68 throw new DimensionMismatchException(yval.length, fval[i].length);
69 }
70
71 final int nI = i + 1;
72 final int pI = i - 1;
73
74 final double nX = xval[nI];
75 final double pX = xval[pI];
76
77 final double deltaX = nX - pX;
78
79 for (int j = 1; j < yLen - 1; j++) {
80 if (zval.length != fval[i][j].length) {
81 throw new DimensionMismatchException(zval.length, fval[i][j].length);
82 }
83
84 final int nJ = j + 1;
85 final int pJ = j - 1;
86
87 final double nY = yval[nJ];
88 final double pY = yval[pJ];
89
90 final double deltaY = nY - pY;
91 final double deltaXY = deltaX * deltaY;
92
93 for (int k = 1; k < zLen - 1; k++) {
94 final int nK = k + 1;
95 final int pK = k - 1;
96
97 final double nZ = zval[nK];
98 final double pZ = zval[pK];
99
100 final double deltaZ = nZ - pZ;
101
102 dFdX[i][j][k] = (fval[nI][j][k] - fval[pI][j][k]) / deltaX;
103 dFdY[i][j][k] = (fval[i][nJ][k] - fval[i][pJ][k]) / deltaY;
104 dFdZ[i][j][k] = (fval[i][j][nK] - fval[i][j][pK]) / deltaZ;
105
106 final double deltaXZ = deltaX * deltaZ;
107 final double deltaYZ = deltaY * deltaZ;
108
109 d2FdXdY[i][j][k] = (fval[nI][nJ][k] - fval[nI][pJ][k] - fval[pI][nJ][k] + fval[pI][pJ][k]) / deltaXY;
110 d2FdXdZ[i][j][k] = (fval[nI][j][nK] - fval[nI][j][pK] - fval[pI][j][nK] + fval[pI][j][pK]) / deltaXZ;
111 d2FdYdZ[i][j][k] = (fval[i][nJ][nK] - fval[i][nJ][pK] - fval[i][pJ][nK] + fval[i][pJ][pK]) / deltaYZ;
112
113 final double deltaXYZ = deltaXY * deltaZ;
114
115 d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] -
116 fval[pI][nJ][nK] + fval[pI][pJ][nK] -
117 fval[nI][nJ][pK] + fval[nI][pJ][pK] +
118 fval[pI][nJ][pK] - fval[pI][pJ][pK]) / deltaXYZ;
119 }
120 }
121 }
122
123
124 return new TricubicInterpolatingFunction(xval, yval, zval, fval,
125 dFdX, dFdY, dFdZ,
126 d2FdXdY, d2FdXdZ, d2FdYdZ,
127 d3FdXdYdZ) {
128
129 @Override
130 public boolean isValidPoint(double x, double y, double z) {
131 return !(x < xval[1] ||
132 x > xval[xval.length - 2] ||
133 y < yval[1] ||
134 y > yval[yval.length - 2] ||
135 z < zval[1] ||
136 z > zval[zval.length - 2]);
137 }
138 };
139 }
140 }