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.math4.legacy.analysis.interpolation; 018 019import org.apache.commons.math4.legacy.exception.DimensionMismatchException; 020import org.apache.commons.math4.legacy.exception.NoDataException; 021import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException; 022import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException; 023import org.apache.commons.math4.legacy.core.MathArrays; 024 025/** 026 * Generates a tricubic interpolating function. 027 * 028 * @since 3.4 029 */ 030public class TricubicInterpolator 031 implements TrivariateGridInterpolator { 032 /** 033 * {@inheritDoc} 034 */ 035 @Override 036 public TricubicInterpolatingFunction interpolate(final double[] xval, 037 final double[] yval, 038 final double[] zval, 039 final double[][][] fval) 040 throws NoDataException, NumberIsTooSmallException, 041 DimensionMismatchException, NonMonotonicSequenceException { 042 if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) { 043 throw new NoDataException(); 044 } 045 if (xval.length != fval.length) { 046 throw new DimensionMismatchException(xval.length, fval.length); 047 } 048 049 MathArrays.checkOrder(xval); 050 MathArrays.checkOrder(yval); 051 MathArrays.checkOrder(zval); 052 053 final int xLen = xval.length; 054 final int yLen = yval.length; 055 final int zLen = zval.length; 056 057 // Approximation to the partial derivatives using finite differences. 058 final double[][][] dFdX = new double[xLen][yLen][zLen]; 059 final double[][][] dFdY = new double[xLen][yLen][zLen]; 060 final double[][][] dFdZ = new double[xLen][yLen][zLen]; 061 final double[][][] d2FdXdY = new double[xLen][yLen][zLen]; 062 final double[][][] d2FdXdZ = new double[xLen][yLen][zLen]; 063 final double[][][] d2FdYdZ = new double[xLen][yLen][zLen]; 064 final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen]; 065 066 for (int i = 1; i < xLen - 1; i++) { 067 if (yval.length != fval[i].length) { 068 throw new DimensionMismatchException(yval.length, fval[i].length); 069 } 070 071 final int nI = i + 1; 072 final int pI = i - 1; 073 074 final double nX = xval[nI]; 075 final double pX = xval[pI]; 076 077 final double deltaX = nX - pX; 078 079 for (int j = 1; j < yLen - 1; j++) { 080 if (zval.length != fval[i][j].length) { 081 throw new DimensionMismatchException(zval.length, fval[i][j].length); 082 } 083 084 final int nJ = j + 1; 085 final int pJ = j - 1; 086 087 final double nY = yval[nJ]; 088 final double pY = yval[pJ]; 089 090 final double deltaY = nY - pY; 091 final double deltaXY = deltaX * deltaY; 092 093 for (int k = 1; k < zLen - 1; k++) { 094 final int nK = k + 1; 095 final int pK = k - 1; 096 097 final double nZ = zval[nK]; 098 final double pZ = zval[pK]; 099 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 // Create the interpolating function. 124 return new TricubicInterpolatingFunction(xval, yval, zval, fval, 125 dFdX, dFdY, dFdZ, 126 d2FdXdY, d2FdXdZ, d2FdYdZ, 127 d3FdXdYdZ) { 128 /** {@inheritDoc} */ 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}