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 * http://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.apache.commons.math.analysis.interpolation; 18 19 import org.apache.commons.math.exception.DimensionMismatchException; 20 import org.apache.commons.math.exception.NoDataException; 21 import org.apache.commons.math.MathException; 22 import org.apache.commons.math.util.MathUtils; 23 24 /** 25 * Generates a tricubic interpolating function. 26 * 27 * @version $Revision$ $Date$ 28 * @since 2.2 29 */ 30 public class TricubicSplineInterpolator 31 implements TrivariateRealGridInterpolator { 32 /** 33 * {@inheritDoc} 34 */ 35 public TricubicSplineInterpolatingFunction interpolate(final double[] xval, 36 final double[] yval, 37 final double[] zval, 38 final double[][][] fval) 39 throws MathException { 40 if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) { 41 throw new NoDataException(); 42 } 43 if (xval.length != fval.length) { 44 throw new DimensionMismatchException(xval.length, fval.length); 45 } 46 47 MathUtils.checkOrder(xval); 48 MathUtils.checkOrder(yval); 49 MathUtils.checkOrder(zval); 50 51 final int xLen = xval.length; 52 final int yLen = yval.length; 53 final int zLen = zval.length; 54 55 // Samples, re-ordered as (z, x, y) and (y, z, x) tuplets 56 // fvalXY[k][i][j] = f(xval[i], yval[j], zval[k]) 57 // fvalZX[j][k][i] = f(xval[i], yval[j], zval[k]) 58 final double[][][] fvalXY = new double[zLen][xLen][yLen]; 59 final double[][][] fvalZX = new double[yLen][zLen][xLen]; 60 for (int i = 0; i < xLen; i++) { 61 if (fval[i].length != yLen) { 62 throw new DimensionMismatchException(fval[i].length, yLen); 63 } 64 65 for (int j = 0; j < yLen; j++) { 66 if (fval[i][j].length != zLen) { 67 throw new DimensionMismatchException(fval[i][j].length, zLen); 68 } 69 70 for (int k = 0; k < zLen; k++) { 71 final double v = fval[i][j][k]; 72 fvalXY[k][i][j] = v; 73 fvalZX[j][k][i] = v; 74 } 75 } 76 } 77 78 final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator(); 79 80 // For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z 81 final BicubicSplineInterpolatingFunction[] xSplineYZ 82 = new BicubicSplineInterpolatingFunction[xLen]; 83 for (int i = 0; i < xLen; i++) { 84 xSplineYZ[i] = bsi.interpolate(yval, zval, fval[i]); 85 } 86 87 // For each line y[j] (0 <= j < yLen), construct a 2D spline in z and x 88 final BicubicSplineInterpolatingFunction[] ySplineZX 89 = new BicubicSplineInterpolatingFunction[yLen]; 90 for (int j = 0; j < yLen; j++) { 91 ySplineZX[j] = bsi.interpolate(zval, xval, fvalZX[j]); 92 } 93 94 // For each line z[k] (0 <= k < zLen), construct a 2D spline in x and y 95 final BicubicSplineInterpolatingFunction[] zSplineXY 96 = new BicubicSplineInterpolatingFunction[zLen]; 97 for (int k = 0; k < zLen; k++) { 98 zSplineXY[k] = bsi.interpolate(xval, yval, fvalXY[k]); 99 } 100 101 // Partial derivatives wrt x and wrt y 102 final double[][][] dFdX = new double[xLen][yLen][zLen]; 103 final double[][][] dFdY = new double[xLen][yLen][zLen]; 104 final double[][][] d2FdXdY = new double[xLen][yLen][zLen]; 105 for (int k = 0; k < zLen; k++) { 106 final BicubicSplineInterpolatingFunction f = zSplineXY[k]; 107 for (int i = 0; i < xLen; i++) { 108 final double x = xval[i]; 109 for (int j = 0; j < yLen; j++) { 110 final double y = yval[j]; 111 dFdX[i][j][k] = f.partialDerivativeX(x, y); 112 dFdY[i][j][k] = f.partialDerivativeY(x, y); 113 d2FdXdY[i][j][k] = f.partialDerivativeXY(x, y); 114 } 115 } 116 } 117 118 // Partial derivatives wrt y and wrt z 119 final double[][][] dFdZ = new double[xLen][yLen][zLen]; 120 final double[][][] d2FdYdZ = new double[xLen][yLen][zLen]; 121 for (int i = 0; i < xLen; i++) { 122 final BicubicSplineInterpolatingFunction f = xSplineYZ[i]; 123 for (int j = 0; j < yLen; j++) { 124 final double y = yval[j]; 125 for (int k = 0; k < zLen; k++) { 126 final double z = zval[k]; 127 dFdZ[i][j][k] = f.partialDerivativeY(y, z); 128 d2FdYdZ[i][j][k] = f.partialDerivativeXY(y, z); 129 } 130 } 131 } 132 133 // Partial derivatives wrt x and wrt z 134 final double[][][] d2FdZdX = new double[xLen][yLen][zLen]; 135 for (int j = 0; j < yLen; j++) { 136 final BicubicSplineInterpolatingFunction f = ySplineZX[j]; 137 for (int k = 0; k < zLen; k++) { 138 final double z = zval[k]; 139 for (int i = 0; i < xLen; i++) { 140 final double x = xval[i]; 141 d2FdZdX[i][j][k] = f.partialDerivativeXY(z, x); 142 } 143 } 144 } 145 146 // Third partial cross-derivatives 147 final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen]; 148 for (int i = 0; i < xLen ; i++) { 149 final int nI = nextIndex(i, xLen); 150 final int pI = previousIndex(i); 151 for (int j = 0; j < yLen; j++) { 152 final int nJ = nextIndex(j, yLen); 153 final int pJ = previousIndex(j); 154 for (int k = 0; k < zLen; k++) { 155 final int nK = nextIndex(k, zLen); 156 final int pK = previousIndex(k); 157 158 // XXX Not sure about this formula 159 d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] - 160 fval[pI][nJ][nK] + fval[pI][pJ][nK] - 161 fval[nI][nJ][pK] + fval[nI][pJ][pK] + 162 fval[pI][nJ][pK] - fval[pI][pJ][pK]) / 163 ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]) * (zval[nK] - zval[pK])) ; 164 } 165 } 166 } 167 168 // Create the interpolating splines 169 return new TricubicSplineInterpolatingFunction(xval, yval, zval, fval, 170 dFdX, dFdY, dFdZ, 171 d2FdXdY, d2FdZdX, d2FdYdZ, 172 d3FdXdYdZ); 173 } 174 175 /** 176 * Compute the next index of an array, clipping if necessary. 177 * It is assumed (but not checked) that {@code i} is larger than or equal to 0}. 178 * 179 * @param i Index 180 * @param max Upper limit of the array 181 * @return the next index 182 */ 183 private int nextIndex(int i, int max) { 184 final int index = i + 1; 185 return index < max ? index : index - 1; 186 } 187 /** 188 * Compute the previous index of an array, clipping if necessary. 189 * It is assumed (but not checked) that {@code i} is smaller than the size of the array. 190 * 191 * @param i Index 192 * @return the previous index 193 */ 194 private int previousIndex(int i) { 195 final int index = i - 1; 196 return index >= 0 ? index : 0; 197 } 198 } 199
