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.analysis.UnivariateRealFunction;
024 import org.apache.commons.math.analysis.BivariateRealFunction;
025 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
026
027 /**
028 * Generates a bicubic interpolation function.
029 * Before interpolating, smoothing of the input data is performed using
030 * splines.
031 * See <b>Handbook on splines for the user</b>, ISBN 084939404X,
032 * chapter 2.
033 *
034 * @version $Revision$ $Date$
035 * @since 2.1
036 */
037 public class SmoothingBicubicSplineInterpolator
038 implements BivariateRealGridInterpolator {
039 /**
040 * {@inheritDoc}
041 */
042 public BivariateRealFunction interpolate(final double[] xval,
043 final double[] yval,
044 final double[][] zval)
045 throws MathException, IllegalArgumentException {
046 if (xval.length == 0 || yval.length == 0 || zval.length == 0) {
047 throw MathRuntimeException.createIllegalArgumentException("no data");
048 }
049 if (xval.length != zval.length) {
050 throw new DimensionMismatchException(xval.length, zval.length);
051 }
052
053 MathUtils.checkOrder(xval, 1, true);
054 MathUtils.checkOrder(yval, 1, true);
055
056 final int xLen = xval.length;
057 final int yLen = yval.length;
058
059 // Samples (first index is y-coordinate, i.e. subarray variable is x)
060 // 0 <= i < xval.length
061 // 0 <= j < yval.length
062 // zX[j][i] = f(xval[i], yval[j])
063 final double[][] zX = new double[yLen][xLen];
064 for (int i = 0; i < xLen; i++) {
065 if (zval[i].length != yLen) {
066 throw new DimensionMismatchException(zval[i].length, yLen);
067 }
068
069 for (int j = 0; j < yLen; j++) {
070 zX[j][i] = zval[i][j];
071 }
072 }
073
074 final SplineInterpolator spInterpolator = new SplineInterpolator();
075
076 // For each line y[j] (0 <= j < yLen), construct a 1D spline with
077 // respect to variable x
078 final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen];
079 for (int j = 0; j < yLen; j++) {
080 ySplineX[j] = spInterpolator.interpolate(xval, zX[j]);
081 }
082
083 // For every knot (xval[i], yval[j]) of the grid, calculate corrected
084 // values zY_1
085 final double[][] zY_1 = new double[xLen][yLen];
086 for (int j = 0; j < yLen; j++) {
087 final PolynomialSplineFunction f = ySplineX[j];
088 for (int i = 0; i < xLen; i++) {
089 zY_1[i][j] = f.value(xval[i]);
090 }
091 }
092
093 // For each line x[i] (0 <= i < xLen), construct a 1D spline with
094 // respect to variable y generated by array zY_1[i]
095 final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen];
096 for (int i = 0; i < xLen; i++) {
097 xSplineY[i] = spInterpolator.interpolate(yval, zY_1[i]);
098 }
099
100 // For every knot (xval[i], yval[j]) of the grid, calculate corrected
101 // values zY_2
102 final double[][] zY_2 = new double[xLen][yLen];
103 for (int i = 0; i < xLen; i++) {
104 final PolynomialSplineFunction f = xSplineY[i];
105 for (int j = 0; j < yLen; j++) {
106 zY_2[i][j] = f.value(yval[j]);
107 }
108 }
109
110 // Partial derivatives with respect to x at the grid knots
111 final double[][] dZdX = new double[xLen][yLen];
112 for (int j = 0; j < yLen; j++) {
113 final UnivariateRealFunction f = ySplineX[j].derivative();
114 for (int i = 0; i < xLen; i++) {
115 dZdX[i][j] = f.value(xval[i]);
116 }
117 }
118
119 // Partial derivatives with respect to y at the grid knots
120 final double[][] dZdY = new double[xLen][yLen];
121 for (int i = 0; i < xLen; i++) {
122 final UnivariateRealFunction f = xSplineY[i].derivative();
123 for (int j = 0; j < yLen; j++) {
124 dZdY[i][j] = f.value(yval[j]);
125 }
126 }
127
128 // Cross partial derivatives
129 final double[][] dZdXdY = new double[xLen][yLen];
130 for (int i = 0; i < xLen ; i++) {
131 final int nI = nextIndex(i, xLen);
132 final int pI = previousIndex(i);
133 for (int j = 0; j < yLen; j++) {
134 final int nJ = nextIndex(j, yLen);
135 final int pJ = previousIndex(j);
136 dZdXdY[i][j] = (zY_2[nI][nJ] - zY_2[nI][pJ] -
137 zY_2[pI][nJ] + zY_2[pI][pJ]) /
138 ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ])) ;
139 }
140 }
141
142 // Create the interpolating splines
143 return new BicubicSplineInterpolatingFunction(xval, yval, zY_2,
144 dZdX, dZdY, dZdXdY);
145 }
146
147 /**
148 * Compute the next index of an array, clipping if necessary.
149 * It is assumed (but not checked) that {@code i} is larger than or equal to 0}.
150 *
151 * @param i Index
152 * @param max Upper limit of the array
153 * @return the next index
154 */
155 private int nextIndex(int i, int max) {
156 final int index = i + 1;
157 return index < max ? index : index - 1;
158 }
159 /**
160 * Compute the previous index of an array, clipping if necessary.
161 * It is assumed (but not checked) that {@code i} is smaller than the size of the array.
162 *
163 * @param i Index
164 * @return the previous index
165 */
166 private int previousIndex(int i) {
167 final int index = i - 1;
168 return index >= 0 ? index : 0;
169 }
170 }