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.MathRuntimeException; 21 import org.apache.commons.math.MathException; 22 import org.apache.commons.math.util.MathUtils; 23 import org.apache.commons.math.util.MathUtils.OrderDirection; 24 import org.apache.commons.math.analysis.BivariateRealFunction; 25 import org.apache.commons.math.analysis.UnivariateRealFunction; 26 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; 27 import org.apache.commons.math.exception.util.LocalizedFormats; 28 29 /** 30 * Generates a bicubic interpolation function. 31 * Before interpolating, smoothing of the input data is performed using 32 * splines. 33 * See <b>Handbook on splines for the user</b>, ISBN 084939404X, 34 * chapter 2. 35 * 36 * @version $Revision: 1059400 $ $Date: 2011-01-15 20:35:27 +0100 (sam. 15 janv. 2011) $ 37 * @since 2.1 38 * @deprecated This class does not perform smoothing; the name is thus misleading. 39 * Please use {@link org.apache.commons.math.analysis.interpolation.BicubicSplineInterpolator} 40 * instead. If smoothing is desired, a tentative implementation is provided in class 41 * {@link org.apache.commons.math.analysis.interpolation.SmoothingPolynomialBicubicSplineInterpolator}. 42 * This class will be removed in math 3.0. 43 */ 44 @Deprecated 45 public class SmoothingBicubicSplineInterpolator 46 implements BivariateRealGridInterpolator { 47 /** 48 * {@inheritDoc} 49 */ interpolate(final double[] xval, final double[] yval, final double[][] zval)50 public BivariateRealFunction interpolate(final double[] xval, 51 final double[] yval, 52 final double[][] zval) 53 throws MathException, IllegalArgumentException { 54 if (xval.length == 0 || yval.length == 0 || zval.length == 0) { 55 throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.NO_DATA); 56 } 57 if (xval.length != zval.length) { 58 throw new DimensionMismatchException(xval.length, zval.length); 59 } 60 61 MathUtils.checkOrder(xval, OrderDirection.INCREASING, true); 62 MathUtils.checkOrder(yval, OrderDirection.INCREASING, true); 63 64 final int xLen = xval.length; 65 final int yLen = yval.length; 66 67 // Samples (first index is y-coordinate, i.e. subarray variable is x) 68 // 0 <= i < xval.length 69 // 0 <= j < yval.length 70 // zX[j][i] = f(xval[i], yval[j]) 71 final double[][] zX = new double[yLen][xLen]; 72 for (int i = 0; i < xLen; i++) { 73 if (zval[i].length != yLen) { 74 throw new DimensionMismatchException(zval[i].length, yLen); 75 } 76 77 for (int j = 0; j < yLen; j++) { 78 zX[j][i] = zval[i][j]; 79 } 80 } 81 82 final SplineInterpolator spInterpolator = new SplineInterpolator(); 83 84 // For each line y[j] (0 <= j < yLen), construct a 1D spline with 85 // respect to variable x 86 final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen]; 87 for (int j = 0; j < yLen; j++) { 88 ySplineX[j] = spInterpolator.interpolate(xval, zX[j]); 89 } 90 91 // For every knot (xval[i], yval[j]) of the grid, calculate corrected 92 // values zY_1 93 final double[][] zY_1 = new double[xLen][yLen]; 94 for (int j = 0; j < yLen; j++) { 95 final PolynomialSplineFunction f = ySplineX[j]; 96 for (int i = 0; i < xLen; i++) { 97 zY_1[i][j] = f.value(xval[i]); 98 } 99 } 100 101 // For each line x[i] (0 <= i < xLen), construct a 1D spline with 102 // respect to variable y generated by array zY_1[i] 103 final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen]; 104 for (int i = 0; i < xLen; i++) { 105 xSplineY[i] = spInterpolator.interpolate(yval, zY_1[i]); 106 } 107 108 // For every knot (xval[i], yval[j]) of the grid, calculate corrected 109 // values zY_2 110 final double[][] zY_2 = new double[xLen][yLen]; 111 for (int i = 0; i < xLen; i++) { 112 final PolynomialSplineFunction f = xSplineY[i]; 113 for (int j = 0; j < yLen; j++) { 114 zY_2[i][j] = f.value(yval[j]); 115 } 116 } 117 118 // Partial derivatives with respect to x at the grid knots 119 final double[][] dZdX = new double[xLen][yLen]; 120 for (int j = 0; j < yLen; j++) { 121 final UnivariateRealFunction f = ySplineX[j].derivative(); 122 for (int i = 0; i < xLen; i++) { 123 dZdX[i][j] = f.value(xval[i]); 124 } 125 } 126 127 // Partial derivatives with respect to y at the grid knots 128 final double[][] dZdY = new double[xLen][yLen]; 129 for (int i = 0; i < xLen; i++) { 130 final UnivariateRealFunction f = xSplineY[i].derivative(); 131 for (int j = 0; j < yLen; j++) { 132 dZdY[i][j] = f.value(yval[j]); 133 } 134 } 135 136 // Cross partial derivatives 137 final double[][] dZdXdY = new double[xLen][yLen]; 138 for (int i = 0; i < xLen ; i++) { 139 final int nI = nextIndex(i, xLen); 140 final int pI = previousIndex(i); 141 for (int j = 0; j < yLen; j++) { 142 final int nJ = nextIndex(j, yLen); 143 final int pJ = previousIndex(j); 144 dZdXdY[i][j] = (zY_2[nI][nJ] - zY_2[nI][pJ] - 145 zY_2[pI][nJ] + zY_2[pI][pJ]) / 146 ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ])); 147 } 148 } 149 150 // Create the interpolating splines 151 return new BicubicSplineInterpolatingFunction(xval, yval, zY_2, 152 dZdX, dZdY, dZdXdY); 153 } 154 155 /** 156 * Compute the next index of an array, clipping if necessary. 157 * It is assumed (but not checked) that {@code i} is larger than or equal to 0}. 158 * 159 * @param i Index 160 * @param max Upper limit of the array 161 * @return the next index 162 */ nextIndex(int i, int max)163 private int nextIndex(int i, int max) { 164 final int index = i + 1; 165 return index < max ? index : index - 1; 166 } 167 /** 168 * Compute the previous index of an array, clipping if necessary. 169 * It is assumed (but not checked) that {@code i} is smaller than the size of the array. 170 * 171 * @param i Index 172 * @return the previous index 173 */ previousIndex(int i)174 private int previousIndex(int i) { 175 final int index = i - 1; 176 return index >= 0 ? index : 0; 177 } 178 } 179