Raymond | dee0849 | 2015-04-02 10:43:13 -0700 | [diff] [blame] | 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 | */ |
| 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 | */ |
| 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 | */ |
| 174 | private int previousIndex(int i) { |
| 175 | final int index = i - 1; |
| 176 | return index >= 0 ? index : 0; |
| 177 | } |
| 178 | } |