001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.math3.analysis.interpolation; 018 019import org.apache.commons.math3.exception.DimensionMismatchException; 020import org.apache.commons.math3.exception.NoDataException; 021import org.apache.commons.math3.exception.NonMonotonicSequenceException; 022import org.apache.commons.math3.exception.NumberIsTooSmallException; 023import org.apache.commons.math3.util.MathArrays; 024 025/** 026 * Generates a tricubic interpolating function. 027 * 028 * @since 2.2 029 * @deprecated To be removed in 4.0 (see MATH-1166). 030 */ 031@Deprecated 032public class TricubicSplineInterpolator 033 implements TrivariateGridInterpolator { 034 /** 035 * {@inheritDoc} 036 */ 037 public TricubicSplineInterpolatingFunction interpolate(final double[] xval, 038 final double[] yval, 039 final double[] zval, 040 final double[][][] fval) 041 throws NoDataException, NumberIsTooSmallException, 042 DimensionMismatchException, NonMonotonicSequenceException { 043 if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) { 044 throw new NoDataException(); 045 } 046 if (xval.length != fval.length) { 047 throw new DimensionMismatchException(xval.length, fval.length); 048 } 049 050 MathArrays.checkOrder(xval); 051 MathArrays.checkOrder(yval); 052 MathArrays.checkOrder(zval); 053 054 final int xLen = xval.length; 055 final int yLen = yval.length; 056 final int zLen = zval.length; 057 058 // Samples, re-ordered as (z, x, y) and (y, z, x) tuplets 059 // fvalXY[k][i][j] = f(xval[i], yval[j], zval[k]) 060 // fvalZX[j][k][i] = f(xval[i], yval[j], zval[k]) 061 final double[][][] fvalXY = new double[zLen][xLen][yLen]; 062 final double[][][] fvalZX = new double[yLen][zLen][xLen]; 063 for (int i = 0; i < xLen; i++) { 064 if (fval[i].length != yLen) { 065 throw new DimensionMismatchException(fval[i].length, yLen); 066 } 067 068 for (int j = 0; j < yLen; j++) { 069 if (fval[i][j].length != zLen) { 070 throw new DimensionMismatchException(fval[i][j].length, zLen); 071 } 072 073 for (int k = 0; k < zLen; k++) { 074 final double v = fval[i][j][k]; 075 fvalXY[k][i][j] = v; 076 fvalZX[j][k][i] = v; 077 } 078 } 079 } 080 081 final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator(true); 082 083 // For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z 084 final BicubicSplineInterpolatingFunction[] xSplineYZ 085 = new BicubicSplineInterpolatingFunction[xLen]; 086 for (int i = 0; i < xLen; i++) { 087 xSplineYZ[i] = bsi.interpolate(yval, zval, fval[i]); 088 } 089 090 // For each line y[j] (0 <= j < yLen), construct a 2D spline in z and x 091 final BicubicSplineInterpolatingFunction[] ySplineZX 092 = new BicubicSplineInterpolatingFunction[yLen]; 093 for (int j = 0; j < yLen; j++) { 094 ySplineZX[j] = bsi.interpolate(zval, xval, fvalZX[j]); 095 } 096 097 // For each line z[k] (0 <= k < zLen), construct a 2D spline in x and y 098 final BicubicSplineInterpolatingFunction[] zSplineXY 099 = new BicubicSplineInterpolatingFunction[zLen]; 100 for (int k = 0; k < zLen; k++) { 101 zSplineXY[k] = bsi.interpolate(xval, yval, fvalXY[k]); 102 } 103 104 // Partial derivatives wrt x and wrt y 105 final double[][][] dFdX = new double[xLen][yLen][zLen]; 106 final double[][][] dFdY = new double[xLen][yLen][zLen]; 107 final double[][][] d2FdXdY = new double[xLen][yLen][zLen]; 108 for (int k = 0; k < zLen; k++) { 109 final BicubicSplineInterpolatingFunction f = zSplineXY[k]; 110 for (int i = 0; i < xLen; i++) { 111 final double x = xval[i]; 112 for (int j = 0; j < yLen; j++) { 113 final double y = yval[j]; 114 dFdX[i][j][k] = f.partialDerivativeX(x, y); 115 dFdY[i][j][k] = f.partialDerivativeY(x, y); 116 d2FdXdY[i][j][k] = f.partialDerivativeXY(x, y); 117 } 118 } 119 } 120 121 // Partial derivatives wrt y and wrt z 122 final double[][][] dFdZ = new double[xLen][yLen][zLen]; 123 final double[][][] d2FdYdZ = new double[xLen][yLen][zLen]; 124 for (int i = 0; i < xLen; i++) { 125 final BicubicSplineInterpolatingFunction f = xSplineYZ[i]; 126 for (int j = 0; j < yLen; j++) { 127 final double y = yval[j]; 128 for (int k = 0; k < zLen; k++) { 129 final double z = zval[k]; 130 dFdZ[i][j][k] = f.partialDerivativeY(y, z); 131 d2FdYdZ[i][j][k] = f.partialDerivativeXY(y, z); 132 } 133 } 134 } 135 136 // Partial derivatives wrt x and wrt z 137 final double[][][] d2FdZdX = new double[xLen][yLen][zLen]; 138 for (int j = 0; j < yLen; j++) { 139 final BicubicSplineInterpolatingFunction f = ySplineZX[j]; 140 for (int k = 0; k < zLen; k++) { 141 final double z = zval[k]; 142 for (int i = 0; i < xLen; i++) { 143 final double x = xval[i]; 144 d2FdZdX[i][j][k] = f.partialDerivativeXY(z, x); 145 } 146 } 147 } 148 149 // Third partial cross-derivatives 150 final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen]; 151 for (int i = 0; i < xLen ; i++) { 152 final int nI = nextIndex(i, xLen); 153 final int pI = previousIndex(i); 154 for (int j = 0; j < yLen; j++) { 155 final int nJ = nextIndex(j, yLen); 156 final int pJ = previousIndex(j); 157 for (int k = 0; k < zLen; k++) { 158 final int nK = nextIndex(k, zLen); 159 final int pK = previousIndex(k); 160 161 // XXX Not sure about this formula 162 d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] - 163 fval[pI][nJ][nK] + fval[pI][pJ][nK] - 164 fval[nI][nJ][pK] + fval[nI][pJ][pK] + 165 fval[pI][nJ][pK] - fval[pI][pJ][pK]) / 166 ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]) * (zval[nK] - zval[pK])) ; 167 } 168 } 169 } 170 171 // Create the interpolating splines 172 return new TricubicSplineInterpolatingFunction(xval, yval, zval, fval, 173 dFdX, dFdY, dFdZ, 174 d2FdXdY, d2FdZdX, d2FdYdZ, 175 d3FdXdYdZ); 176 } 177 178 /** 179 * Compute the next index of an array, clipping if necessary. 180 * It is assumed (but not checked) that {@code i} is larger than or equal to 0. 181 * 182 * @param i Index 183 * @param max Upper limit of the array 184 * @return the next index 185 */ 186 private int nextIndex(int i, int max) { 187 final int index = i + 1; 188 return index < max ? index : index - 1; 189 } 190 /** 191 * Compute the previous index of an array, clipping if necessary. 192 * It is assumed (but not checked) that {@code i} is smaller than the size of the array. 193 * 194 * @param i Index 195 * @return the previous index 196 */ 197 private int previousIndex(int i) { 198 final int index = i - 1; 199 return index >= 0 ? index : 0; 200 } 201}