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 */ 017 package org.apache.commons.math.analysis.interpolation; 018 019 import org.apache.commons.math.DimensionMismatchException; 020 import org.apache.commons.math.MathRuntimeException; 021 import org.apache.commons.math.MathException; 022 import org.apache.commons.math.util.MathUtils; 023 import org.apache.commons.math.util.MathUtils.OrderDirection; 024 import org.apache.commons.math.analysis.BivariateRealFunction; 025 import org.apache.commons.math.analysis.UnivariateRealFunction; 026 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; 027 import org.apache.commons.math.exception.util.LocalizedFormats; 028 029 /** 030 * Generates a bicubic interpolation function. 031 * Before interpolating, smoothing of the input data is performed using 032 * splines. 033 * See <b>Handbook on splines for the user</b>, ISBN 084939404X, 034 * chapter 2. 035 * 036 * @version $Revision: 1059400 $ $Date: 2011-01-15 20:35:27 +0100 (sam. 15 janv. 2011) $ 037 * @since 2.1 038 * @deprecated This class does not perform smoothing; the name is thus misleading. 039 * Please use {@link org.apache.commons.math.analysis.interpolation.BicubicSplineInterpolator} 040 * instead. If smoothing is desired, a tentative implementation is provided in class 041 * {@link org.apache.commons.math.analysis.interpolation.SmoothingPolynomialBicubicSplineInterpolator}. 042 * This class will be removed in math 3.0. 043 */ 044 @Deprecated 045 public class SmoothingBicubicSplineInterpolator 046 implements BivariateRealGridInterpolator { 047 /** 048 * {@inheritDoc} 049 */ 050 public BivariateRealFunction interpolate(final double[] xval, 051 final double[] yval, 052 final double[][] zval) 053 throws MathException, IllegalArgumentException { 054 if (xval.length == 0 || yval.length == 0 || zval.length == 0) { 055 throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.NO_DATA); 056 } 057 if (xval.length != zval.length) { 058 throw new DimensionMismatchException(xval.length, zval.length); 059 } 060 061 MathUtils.checkOrder(xval, OrderDirection.INCREASING, true); 062 MathUtils.checkOrder(yval, OrderDirection.INCREASING, true); 063 064 final int xLen = xval.length; 065 final int yLen = yval.length; 066 067 // Samples (first index is y-coordinate, i.e. subarray variable is x) 068 // 0 <= i < xval.length 069 // 0 <= j < yval.length 070 // zX[j][i] = f(xval[i], yval[j]) 071 final double[][] zX = new double[yLen][xLen]; 072 for (int i = 0; i < xLen; i++) { 073 if (zval[i].length != yLen) { 074 throw new DimensionMismatchException(zval[i].length, yLen); 075 } 076 077 for (int j = 0; j < yLen; j++) { 078 zX[j][i] = zval[i][j]; 079 } 080 } 081 082 final SplineInterpolator spInterpolator = new SplineInterpolator(); 083 084 // For each line y[j] (0 <= j < yLen), construct a 1D spline with 085 // respect to variable x 086 final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen]; 087 for (int j = 0; j < yLen; j++) { 088 ySplineX[j] = spInterpolator.interpolate(xval, zX[j]); 089 } 090 091 // For every knot (xval[i], yval[j]) of the grid, calculate corrected 092 // values zY_1 093 final double[][] zY_1 = new double[xLen][yLen]; 094 for (int j = 0; j < yLen; j++) { 095 final PolynomialSplineFunction f = ySplineX[j]; 096 for (int i = 0; i < xLen; i++) { 097 zY_1[i][j] = f.value(xval[i]); 098 } 099 } 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 }