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.util.MathUtils; 020 import org.apache.commons.math.MathRuntimeException; 021 import org.apache.commons.math.DimensionMismatchException; 022 import org.apache.commons.math.analysis.BivariateRealFunction; 023 024 /** 025 * Function that implements the 026 * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation"> 027 * bicubic spline interpolation</a>. 028 * 029 * @version $Revision$ $Date$ 030 * @since 2.1 031 */ 032 public class BicubicSplineInterpolatingFunction 033 implements BivariateRealFunction { 034 /** 035 * Matrix to compute the spline coefficients from the function values 036 * and function derivatives values 037 */ 038 private static final double[][] AINV = { 039 { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, 040 { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 }, 041 { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 }, 042 { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 }, 043 { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 }, 044 { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 }, 045 { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 }, 046 { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 }, 047 { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 }, 048 { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 }, 049 { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 }, 050 { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 }, 051 { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 }, 052 { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 }, 053 { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 }, 054 { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 } 055 }; 056 057 /** Samples x-coordinates */ 058 private final double[] xval; 059 /** Samples y-coordinates */ 060 private final double[] yval; 061 /** Set of cubic splines pacthing the whole data grid */ 062 private final BicubicSplineFunction[][] splines; 063 064 /** 065 * @param x Sample values of the x-coordinate, in increasing order 066 * @param y Sample values of the y-coordinate, in increasing order 067 * @param z Values of the function on every grid point 068 * @param dZdX Values of the partial derivative of function with respect 069 * to x on every grid point 070 * @param dZdY Values of the partial derivative of function with respect 071 * to y on every grid point 072 * @param dZdXdY Values of the cross partial derivative of function on 073 * every grid point 074 * @throws DimensionMismatchException if the various arrays do not contain 075 * the expected number of elements. 076 * @throws IllegalArgumentException if {@code x} or {@code y} are not strictly 077 * increasing. 078 */ 079 public BicubicSplineInterpolatingFunction(double[] x, 080 double[] y, 081 double[][] z, 082 double[][] dZdX, 083 double[][] dZdY, 084 double[][] dZdXdY) 085 throws DimensionMismatchException { 086 final int xLen = x.length; 087 final int yLen = y.length; 088 089 if (xLen == 0 || yLen == 0 || z.length == 0 || z[0].length == 0) { 090 throw MathRuntimeException.createIllegalArgumentException("no data"); 091 } 092 if (xLen != z.length) { 093 throw new DimensionMismatchException(xLen, z.length); 094 } 095 if (xLen != dZdX.length) { 096 throw new DimensionMismatchException(xLen, dZdX.length); 097 } 098 if (xLen != dZdY.length) { 099 throw new DimensionMismatchException(xLen, dZdY.length); 100 } 101 if (xLen != dZdXdY.length) { 102 throw new DimensionMismatchException(xLen, dZdXdY.length); 103 } 104 105 MathUtils.checkOrder(x, 1, true); 106 MathUtils.checkOrder(y, 1, true); 107 108 xval = x.clone(); 109 yval = y.clone(); 110 111 final int lastI = xLen - 1; 112 final int lastJ = yLen - 1; 113 splines = new BicubicSplineFunction[lastI][lastJ]; 114 115 for (int i = 0; i < lastI; i++) { 116 if (z[i].length != yLen) { 117 throw new DimensionMismatchException(z[i].length, yLen); 118 } 119 if (dZdX[i].length != yLen) { 120 throw new DimensionMismatchException(dZdX[i].length, yLen); 121 } 122 if (dZdY[i].length != yLen) { 123 throw new DimensionMismatchException(dZdY[i].length, yLen); 124 } 125 if (dZdXdY[i].length != yLen) { 126 throw new DimensionMismatchException(dZdXdY[i].length, yLen); 127 } 128 final int ip1 = i + 1; 129 for (int j = 0; j < lastJ; j++) { 130 final int jp1 = j + 1; 131 final double[] beta = new double[] { 132 z[i][j], z[ip1][j], z[i][jp1], z[ip1][jp1], 133 dZdX[i][j], dZdX[ip1][j], dZdX[i][jp1], dZdX[ip1][jp1], 134 dZdY[i][j], dZdY[ip1][j], dZdY[i][jp1], dZdY[ip1][jp1], 135 dZdXdY[i][j], dZdXdY[ip1][j], dZdXdY[i][jp1], dZdXdY[ip1][jp1] 136 }; 137 138 splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta)); 139 } 140 } 141 } 142 143 /** 144 * {@inheritDoc} 145 */ 146 public double value(double x, double y) { 147 final int i = searchIndex(x, xval); 148 if (i == -1) { 149 throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range", 150 x, xval[0], xval[xval.length - 1]); 151 } 152 final int j = searchIndex(y, yval); 153 if (j == -1) { 154 throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range", 155 y, yval[0], yval[yval.length - 1]); 156 } 157 158 final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); 159 final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); 160 161 return splines[i][j].value(xN, yN); 162 } 163 164 /** 165 * @param c coordinate 166 * @param val coordinate samples 167 * @return the index in {@code val} corresponding to the interval 168 * containing {@code c}, or {@code -1} if {@code c} is out of the 169 * range defined by the end values of {@code val} 170 */ 171 private int searchIndex(double c, double[] val) { 172 if (c < val[0]) { 173 return -1; 174 } 175 176 int max = val.length; 177 for (int i = 1; i < max; i++) { 178 if (c <= val[i]) { 179 return i - 1; 180 } 181 } 182 183 return -1; 184 } 185 186 /** 187 * Compute the spline coefficients from the list of function values and 188 * function partial derivatives values at the four corners of a grid 189 * element. They must be specified in the following order: 190 * <ul> 191 * <li>f(0,0)</li> 192 * <li>f(1,0)</li> 193 * <li>f(0,1)</li> 194 * <li>f(1,1)</li> 195 * <li>fx(0,0)</li> 196 * <li>fx(1,0)</li> 197 * <li>fx(0,1)</li> 198 * <li>fx(1,1)</li> 199 * <li>fy(0,0)</li> 200 * <li>fy(1,0)</li> 201 * <li>fy(0,1)</li> 202 * <li>fy(1,1)</li> 203 * <li>fxy(0,0)</li> 204 * <li>fxy(1,0)</li> 205 * <li>fxy(0,1)</li> 206 * <li>fxy(1,1)</li> 207 * </ul> 208 * @param beta List of function values and function partial derivatives 209 * values 210 * @return the spline coefficients 211 */ 212 private double[] computeSplineCoefficients(double[] beta) { 213 final double[] a = new double[16]; 214 215 for (int i = 0; i < 16; i++) { 216 double result = 0; 217 final double[] row = AINV[i]; 218 for (int j = 0; j < 16; j++) { 219 result += row[j] * beta[j]; 220 } 221 a[i] = result; 222 } 223 224 return a; 225 } 226 } 227 /** 228 * 2D-spline function. 229 * 230 * @version $Revision$ $Date$ 231 */ 232 class BicubicSplineFunction 233 implements BivariateRealFunction { 234 //CHECKSTYLE: stop MultipleVariableDeclarations 235 /** Coefficients */ 236 private final double 237 a00, a01, a02, a03, 238 a10, a11, a12, a13, 239 a20, a21, a22, a23, 240 a30, a31, a32, a33; 241 //CHECKSTYLE: resume MultipleVariableDeclarations 242 243 /** 244 * @param a Spline coefficients 245 */ 246 public BicubicSplineFunction(double[] a) { 247 this.a00 = a[0]; 248 this.a10 = a[1]; 249 this.a20 = a[2]; 250 this.a30 = a[3]; 251 this.a01 = a[4]; 252 this.a11 = a[5]; 253 this.a21 = a[6]; 254 this.a31 = a[7]; 255 this.a02 = a[8]; 256 this.a12 = a[9]; 257 this.a22 = a[10]; 258 this.a32 = a[11]; 259 this.a03 = a[12]; 260 this.a13 = a[13]; 261 this.a23 = a[14]; 262 this.a33 = a[15]; 263 } 264 265 /** 266 * @param x x-coordinate of the interpolation point 267 * @param y y-coordinate of the interpolation point 268 * @return the interpolated value. 269 */ 270 public double value(double x, double y) { 271 if (x < 0 || x > 1) { 272 throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range", 273 x, 0, 1); 274 } 275 if (y < 0 || y > 1) { 276 throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range", 277 y, 0, 1); 278 } 279 280 final double x2 = x * x; 281 final double x3 = x2 * x; 282 final double y2 = y * y; 283 final double y3 = y2 * y; 284 285 return a00 + a01 * y + a02 * y2 + a03 * y3 + 286 a10 * x + a11 * x * y + a12 * x * y2 + a13 * x * y3 + 287 a20 * x2 + a21 * x2 * y + a22 * x2 * y2 + a23 * x2 * y3 + 288 a30 * x3 + a31 * x3 * y + a32 * x3 * y2 + a33 * x3 * y3; 289 } 290 }