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 */ 017package org.apache.commons.math3.analysis.interpolation; 018 019import java.util.Arrays; 020import org.apache.commons.math3.analysis.BivariateFunction; 021import org.apache.commons.math3.exception.DimensionMismatchException; 022import org.apache.commons.math3.exception.NoDataException; 023import org.apache.commons.math3.exception.OutOfRangeException; 024import org.apache.commons.math3.exception.NonMonotonicSequenceException; 025import org.apache.commons.math3.util.MathArrays; 026 027/** 028 * Function that implements the 029 * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation"> 030 * bicubic spline interpolation</a>. 031 * 032 * @since 3.4 033 */ 034public class BicubicInterpolatingFunction 035 implements BivariateFunction { 036 /** Number of coefficients. */ 037 private static final int NUM_COEFF = 16; 038 /** 039 * Matrix to compute the spline coefficients from the function values 040 * and function derivatives values 041 */ 042 private static final double[][] AINV = { 043 { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, 044 { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 }, 045 { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 }, 046 { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 }, 047 { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 }, 048 { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 }, 049 { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 }, 050 { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 }, 051 { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 }, 052 { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 }, 053 { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 }, 054 { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 }, 055 { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 }, 056 { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 }, 057 { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 }, 058 { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 } 059 }; 060 061 /** Samples x-coordinates */ 062 private final double[] xval; 063 /** Samples y-coordinates */ 064 private final double[] yval; 065 /** Set of cubic splines patching the whole data grid */ 066 private final BicubicFunction[][] splines; 067 068 /** 069 * @param x Sample values of the x-coordinate, in increasing order. 070 * @param y Sample values of the y-coordinate, in increasing order. 071 * @param f Values of the function on every grid point. 072 * @param dFdX Values of the partial derivative of function with respect 073 * to x on every grid point. 074 * @param dFdY Values of the partial derivative of function with respect 075 * to y on every grid point. 076 * @param d2FdXdY Values of the cross partial derivative of function on 077 * every grid point. 078 * @throws DimensionMismatchException if the various arrays do not contain 079 * the expected number of elements. 080 * @throws NonMonotonicSequenceException if {@code x} or {@code y} are 081 * not strictly increasing. 082 * @throws NoDataException if any of the arrays has zero length. 083 */ 084 public BicubicInterpolatingFunction(double[] x, 085 double[] y, 086 double[][] f, 087 double[][] dFdX, 088 double[][] dFdY, 089 double[][] d2FdXdY) 090 throws DimensionMismatchException, 091 NoDataException, 092 NonMonotonicSequenceException { 093 final int xLen = x.length; 094 final int yLen = y.length; 095 096 if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) { 097 throw new NoDataException(); 098 } 099 if (xLen != f.length) { 100 throw new DimensionMismatchException(xLen, f.length); 101 } 102 if (xLen != dFdX.length) { 103 throw new DimensionMismatchException(xLen, dFdX.length); 104 } 105 if (xLen != dFdY.length) { 106 throw new DimensionMismatchException(xLen, dFdY.length); 107 } 108 if (xLen != d2FdXdY.length) { 109 throw new DimensionMismatchException(xLen, d2FdXdY.length); 110 } 111 112 MathArrays.checkOrder(x); 113 MathArrays.checkOrder(y); 114 115 xval = x.clone(); 116 yval = y.clone(); 117 118 final int lastI = xLen - 1; 119 final int lastJ = yLen - 1; 120 splines = new BicubicFunction[lastI][lastJ]; 121 122 for (int i = 0; i < lastI; i++) { 123 if (f[i].length != yLen) { 124 throw new DimensionMismatchException(f[i].length, yLen); 125 } 126 if (dFdX[i].length != yLen) { 127 throw new DimensionMismatchException(dFdX[i].length, yLen); 128 } 129 if (dFdY[i].length != yLen) { 130 throw new DimensionMismatchException(dFdY[i].length, yLen); 131 } 132 if (d2FdXdY[i].length != yLen) { 133 throw new DimensionMismatchException(d2FdXdY[i].length, yLen); 134 } 135 final int ip1 = i + 1; 136 final double xR = xval[ip1] - xval[i]; 137 for (int j = 0; j < lastJ; j++) { 138 final int jp1 = j + 1; 139 final double yR = yval[jp1] - yval[j]; 140 final double xRyR = xR * yR; 141 final double[] beta = new double[] { 142 f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1], 143 dFdX[i][j] * xR, dFdX[ip1][j] * xR, dFdX[i][jp1] * xR, dFdX[ip1][jp1] * xR, 144 dFdY[i][j] * yR, dFdY[ip1][j] * yR, dFdY[i][jp1] * yR, dFdY[ip1][jp1] * yR, 145 d2FdXdY[i][j] * xRyR, d2FdXdY[ip1][j] * xRyR, d2FdXdY[i][jp1] * xRyR, d2FdXdY[ip1][jp1] * xRyR 146 }; 147 148 splines[i][j] = new BicubicFunction(computeSplineCoefficients(beta)); 149 } 150 } 151 } 152 153 /** 154 * {@inheritDoc} 155 */ 156 public double value(double x, double y) 157 throws OutOfRangeException { 158 final int i = searchIndex(x, xval); 159 final int j = searchIndex(y, yval); 160 161 final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); 162 final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); 163 164 return splines[i][j].value(xN, yN); 165 } 166 167 /** 168 * Indicates whether a point is within the interpolation range. 169 * 170 * @param x First coordinate. 171 * @param y Second coordinate. 172 * @return {@code true} if (x, y) is a valid point. 173 */ 174 public boolean isValidPoint(double x, double y) { 175 if (x < xval[0] || 176 x > xval[xval.length - 1] || 177 y < yval[0] || 178 y > yval[yval.length - 1]) { 179 return false; 180 } else { 181 return true; 182 } 183 } 184 185 /** 186 * @param c Coordinate. 187 * @param val Coordinate samples. 188 * @return the index in {@code val} corresponding to the interval 189 * containing {@code c}. 190 * @throws OutOfRangeException if {@code c} is out of the 191 * range defined by the boundary values of {@code val}. 192 */ 193 private int searchIndex(double c, double[] val) { 194 final int r = Arrays.binarySearch(val, c); 195 196 if (r == -1 || 197 r == -val.length - 1) { 198 throw new OutOfRangeException(c, val[0], val[val.length - 1]); 199 } 200 201 if (r < 0) { 202 // "c" in within an interpolation sub-interval: Return the 203 // index of the sample at the lower end of the sub-interval. 204 return -r - 2; 205 } 206 final int last = val.length - 1; 207 if (r == last) { 208 // "c" is the last sample of the range: Return the index 209 // of the sample at the lower end of the last sub-interval. 210 return last - 1; 211 } 212 213 // "c" is another sample point. 214 return r; 215 } 216 217 /** 218 * Compute the spline coefficients from the list of function values and 219 * function partial derivatives values at the four corners of a grid 220 * element. They must be specified in the following order: 221 * <ul> 222 * <li>f(0,0)</li> 223 * <li>f(1,0)</li> 224 * <li>f(0,1)</li> 225 * <li>f(1,1)</li> 226 * <li>f<sub>x</sub>(0,0)</li> 227 * <li>f<sub>x</sub>(1,0)</li> 228 * <li>f<sub>x</sub>(0,1)</li> 229 * <li>f<sub>x</sub>(1,1)</li> 230 * <li>f<sub>y</sub>(0,0)</li> 231 * <li>f<sub>y</sub>(1,0)</li> 232 * <li>f<sub>y</sub>(0,1)</li> 233 * <li>f<sub>y</sub>(1,1)</li> 234 * <li>f<sub>xy</sub>(0,0)</li> 235 * <li>f<sub>xy</sub>(1,0)</li> 236 * <li>f<sub>xy</sub>(0,1)</li> 237 * <li>f<sub>xy</sub>(1,1)</li> 238 * </ul> 239 * where the subscripts indicate the partial derivative with respect to 240 * the corresponding variable(s). 241 * 242 * @param beta List of function values and function partial derivatives 243 * values. 244 * @return the spline coefficients. 245 */ 246 private double[] computeSplineCoefficients(double[] beta) { 247 final double[] a = new double[NUM_COEFF]; 248 249 for (int i = 0; i < NUM_COEFF; i++) { 250 double result = 0; 251 final double[] row = AINV[i]; 252 for (int j = 0; j < NUM_COEFF; j++) { 253 result += row[j] * beta[j]; 254 } 255 a[i] = result; 256 } 257 258 return a; 259 } 260} 261 262/** 263 * Bicubic function. 264 */ 265class BicubicFunction implements BivariateFunction { 266 /** Number of points. */ 267 private static final short N = 4; 268 /** Coefficients */ 269 private final double[][] a; 270 271 /** 272 * Simple constructor. 273 * 274 * @param coeff Spline coefficients. 275 */ 276 BicubicFunction(double[] coeff) { 277 a = new double[N][N]; 278 for (int j = 0; j < N; j++) { 279 final double[] aJ = a[j]; 280 for (int i = 0; i < N; i++) { 281 aJ[i] = coeff[i * N + j]; 282 } 283 } 284 } 285 286 /** 287 * {@inheritDoc} 288 */ 289 public double value(double x, double y) { 290 if (x < 0 || x > 1) { 291 throw new OutOfRangeException(x, 0, 1); 292 } 293 if (y < 0 || y > 1) { 294 throw new OutOfRangeException(y, 0, 1); 295 } 296 297 final double x2 = x * x; 298 final double x3 = x2 * x; 299 final double[] pX = {1, x, x2, x3}; 300 301 final double y2 = y * y; 302 final double y3 = y2 * y; 303 final double[] pY = {1, y, y2, y3}; 304 305 return apply(pX, pY, a); 306 } 307 308 /** 309 * Compute the value of the bicubic polynomial. 310 * 311 * @param pX Powers of the x-coordinate. 312 * @param pY Powers of the y-coordinate. 313 * @param coeff Spline coefficients. 314 * @return the interpolated value. 315 */ 316 private double apply(double[] pX, double[] pY, double[][] coeff) { 317 double result = 0; 318 for (int i = 0; i < N; i++) { 319 final double r = MathArrays.linearCombination(coeff[i], pY); 320 result += r * pX[i]; 321 } 322 323 return result; 324 } 325}