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.analysis.polynomials.PolynomialSplineFunction; 022import org.apache.commons.math3.exception.DimensionMismatchException; 023import org.apache.commons.math3.exception.InsufficientDataException; 024import org.apache.commons.math3.exception.NoDataException; 025import org.apache.commons.math3.exception.NullArgumentException; 026import org.apache.commons.math3.exception.OutOfRangeException; 027import org.apache.commons.math3.exception.NonMonotonicSequenceException; 028import org.apache.commons.math3.util.MathArrays; 029 030/** 031 * Function that implements the 032 * <a href="http://www.paulinternet.nl/?page=bicubic">bicubic spline</a> 033 * interpolation. 034 * This implementation currently uses {@link AkimaSplineInterpolator} as the 035 * underlying one-dimensional interpolator, which requires 5 sample points; 036 * insufficient data will raise an exception when the 037 * {@link #value(double,double) value} method is called. 038 * 039 * @since 3.4 040 */ 041public class PiecewiseBicubicSplineInterpolatingFunction 042 implements BivariateFunction { 043 /** The minimum number of points that are needed to compute the function. */ 044 private static final int MIN_NUM_POINTS = 5; 045 /** Samples x-coordinates */ 046 private final double[] xval; 047 /** Samples y-coordinates */ 048 private final double[] yval; 049 /** Set of cubic splines patching the whole data grid */ 050 private final double[][] fval; 051 052 /** 053 * @param x Sample values of the x-coordinate, in increasing order. 054 * @param y Sample values of the y-coordinate, in increasing order. 055 * @param f Values of the function on every grid point. the expected number 056 * of elements. 057 * @throws NonMonotonicSequenceException if {@code x} or {@code y} are not 058 * strictly increasing. 059 * @throws NullArgumentException if any of the arguments are null 060 * @throws NoDataException if any of the arrays has zero length. 061 * @throws DimensionMismatchException if the length of x and y don't match the row, column 062 * height of f 063 */ 064 public PiecewiseBicubicSplineInterpolatingFunction(double[] x, 065 double[] y, 066 double[][] f) 067 throws DimensionMismatchException, 068 NullArgumentException, 069 NoDataException, 070 NonMonotonicSequenceException { 071 if (x == null || 072 y == null || 073 f == null || 074 f[0] == null) { 075 throw new NullArgumentException(); 076 } 077 078 final int xLen = x.length; 079 final int yLen = y.length; 080 081 if (xLen == 0 || 082 yLen == 0 || 083 f.length == 0 || 084 f[0].length == 0) { 085 throw new NoDataException(); 086 } 087 088 if (xLen < MIN_NUM_POINTS || 089 yLen < MIN_NUM_POINTS || 090 f.length < MIN_NUM_POINTS || 091 f[0].length < MIN_NUM_POINTS) { 092 throw new InsufficientDataException(); 093 } 094 095 if (xLen != f.length) { 096 throw new DimensionMismatchException(xLen, f.length); 097 } 098 099 if (yLen != f[0].length) { 100 throw new DimensionMismatchException(yLen, f[0].length); 101 } 102 103 MathArrays.checkOrder(x); 104 MathArrays.checkOrder(y); 105 106 xval = x.clone(); 107 yval = y.clone(); 108 fval = f.clone(); 109 } 110 111 /** 112 * {@inheritDoc} 113 */ 114 public double value(double x, 115 double y) 116 throws OutOfRangeException { 117 final AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator(); 118 final int offset = 2; 119 final int count = offset + 3; 120 final int i = searchIndex(x, xval, offset, count); 121 final int j = searchIndex(y, yval, offset, count); 122 123 final double xArray[] = new double[count]; 124 final double yArray[] = new double[count]; 125 final double zArray[] = new double[count]; 126 final double interpArray[] = new double[count]; 127 128 for (int index = 0; index < count; index++) { 129 xArray[index] = xval[i + index]; 130 yArray[index] = yval[j + index]; 131 } 132 133 for (int zIndex = 0; zIndex < count; zIndex++) { 134 for (int index = 0; index < count; index++) { 135 zArray[index] = fval[i + index][j + zIndex]; 136 } 137 final PolynomialSplineFunction spline = interpolator.interpolate(xArray, zArray); 138 interpArray[zIndex] = spline.value(x); 139 } 140 141 final PolynomialSplineFunction spline = interpolator.interpolate(yArray, interpArray); 142 143 double returnValue = spline.value(y); 144 145 return returnValue; 146 } 147 148 /** 149 * Indicates whether a point is within the interpolation range. 150 * 151 * @param x First coordinate. 152 * @param y Second coordinate. 153 * @return {@code true} if (x, y) is a valid point. 154 * @since 3.3 155 */ 156 public boolean isValidPoint(double x, 157 double y) { 158 if (x < xval[0] || 159 x > xval[xval.length - 1] || 160 y < yval[0] || 161 y > yval[yval.length - 1]) { 162 return false; 163 } else { 164 return true; 165 } 166 } 167 168 /** 169 * @param c Coordinate. 170 * @param val Coordinate samples. 171 * @param offset how far back from found value to offset for querying 172 * @param count total number of elements forward from beginning that will be 173 * queried 174 * @return the index in {@code val} corresponding to the interval containing 175 * {@code c}. 176 * @throws OutOfRangeException if {@code c} is out of the range defined by 177 * the boundary values of {@code val}. 178 */ 179 private int searchIndex(double c, 180 double[] val, 181 int offset, 182 int count) { 183 int r = Arrays.binarySearch(val, c); 184 185 if (r == -1 || r == -val.length - 1) { 186 throw new OutOfRangeException(c, val[0], val[val.length - 1]); 187 } 188 189 if (r < 0) { 190 // "c" in within an interpolation sub-interval, which returns 191 // negative 192 // need to remove the negative sign for consistency 193 r = -r - offset - 1; 194 } else { 195 r -= offset; 196 } 197 198 if (r < 0) { 199 r = 0; 200 } 201 202 if ((r + count) >= val.length) { 203 // "c" is the last sample of the range: Return the index 204 // of the sample at the lower end of the last sub-interval. 205 r = val.length - count; 206 } 207 208 return r; 209 } 210}