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.DimensionMismatchException; 20 import org.apache.commons.math.MathException; 21 import org.apache.commons.math.analysis.UnivariateRealFunction; 22 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; 23 import org.apache.commons.math.exception.NoDataException; 24 import org.apache.commons.math.util.MathUtils; 25 26 /** 27 * Generates a bicubic interpolating function. 28 * 29 * @version $Revision: 980944 $ $Date: 2010-07-30 22:31:11 +0200 (ven. 30 juil. 2010) $ 30 * @since 2.2 31 */ 32 public class BicubicSplineInterpolator 33 implements BivariateRealGridInterpolator { 34 /** 35 * {@inheritDoc} 36 */ 37 public BicubicSplineInterpolatingFunction interpolate(final double[] xval, 38 final double[] yval, 39 final double[][] fval) 40 throws MathException, IllegalArgumentException { 41 if (xval.length == 0 || yval.length == 0 || fval.length == 0) { 42 throw new NoDataException(); 43 } 44 if (xval.length != fval.length) { 45 throw new DimensionMismatchException(xval.length, fval.length); 46 } 47 48 MathUtils.checkOrder(xval); 49 MathUtils.checkOrder(yval); 50 51 final int xLen = xval.length; 52 final int yLen = yval.length; 53 54 // Samples (first index is y-coordinate, i.e. subarray variable is x) 55 // 0 <= i < xval.length 56 // 0 <= j < yval.length 57 // fX[j][i] = f(xval[i], yval[j]) 58 final double[][] fX = new double[yLen][xLen]; 59 for (int i = 0; i < xLen; i++) { 60 if (fval[i].length != yLen) { 61 throw new DimensionMismatchException(fval[i].length, yLen); 62 } 63 64 for (int j = 0; j < yLen; j++) { 65 fX[j][i] = fval[i][j]; 66 } 67 } 68 69 final SplineInterpolator spInterpolator = new SplineInterpolator(); 70 71 // For each line y[j] (0 <= j < yLen), construct a 1D spline with 72 // respect to variable x 73 final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen]; 74 for (int j = 0; j < yLen; j++) { 75 ySplineX[j] = spInterpolator.interpolate(xval, fX[j]); 76 } 77 78 // For each line x[i] (0 <= i < xLen), construct a 1D spline with 79 // respect to variable y generated by array fY_1[i] 80 final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen]; 81 for (int i = 0; i < xLen; i++) { 82 xSplineY[i] = spInterpolator.interpolate(yval, fval[i]); 83 } 84 85 // Partial derivatives with respect to x at the grid knots 86 final double[][] dFdX = new double[xLen][yLen]; 87 for (int j = 0; j < yLen; j++) { 88 final UnivariateRealFunction f = ySplineX[j].derivative(); 89 for (int i = 0; i < xLen; i++) { 90 dFdX[i][j] = f.value(xval[i]); 91 } 92 } 93 94 // Partial derivatives with respect to y at the grid knots 95 final double[][] dFdY = new double[xLen][yLen]; 96 for (int i = 0; i < xLen; i++) { 97 final UnivariateRealFunction f = xSplineY[i].derivative(); 98 for (int j = 0; j < yLen; j++) { 99 dFdY[i][j] = f.value(yval[j]); 100 } 101 } 102 103 // Cross partial derivatives 104 final double[][] d2FdXdY = new double[xLen][yLen]; 105 for (int i = 0; i < xLen ; i++) { 106 final int nI = nextIndex(i, xLen); 107 final int pI = previousIndex(i); 108 for (int j = 0; j < yLen; j++) { 109 final int nJ = nextIndex(j, yLen); 110 final int pJ = previousIndex(j); 111 d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] - 112 fval[pI][nJ] + fval[pI][pJ]) / 113 ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ])); 114 } 115 } 116 117 // Create the interpolating splines 118 return new BicubicSplineInterpolatingFunction(xval, yval, fval, 119 dFdX, dFdY, d2FdXdY); 120 } 121 122 /** 123 * Compute the next index of an array, clipping if necessary. 124 * It is assumed (but not checked) that {@code i} is larger than or equal to 0}. 125 * 126 * @param i Index 127 * @param max Upper limit of the array 128 * @return the next index 129 */ 130 private int nextIndex(int i, int max) { 131 final int index = i + 1; 132 return index < max ? index : index - 1; 133 } 134 /** 135 * Compute the previous index of an array, clipping if necessary. 136 * It is assumed (but not checked) that {@code i} is smaller than the size of the array. 137 * 138 * @param i Index 139 * @return the previous index 140 */ 141 private int previousIndex(int i) { 142 final int index = i - 1; 143 return index >= 0 ? index : 0; 144 } 145 } 146
