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.exception.DimensionMismatchException;
020 import org.apache.commons.math.exception.NoDataException;
021 import org.apache.commons.math.MathException;
022 import org.apache.commons.math.util.MathUtils;
023 import org.apache.commons.math.optimization.general.GaussNewtonOptimizer;
024 import org.apache.commons.math.optimization.fitting.PolynomialFitter;
025 import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
026
027 /**
028 * Generates a bicubic interpolation function.
029 * Prior to generating the interpolating function, the input is smoothed using
030 * polynomial fitting.
031 *
032 * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $
033 * @since 2.2
034 */
035 public class SmoothingPolynomialBicubicSplineInterpolator
036 extends BicubicSplineInterpolator {
037
038 /** Fitter for x. */
039 private final PolynomialFitter xFitter;
040
041 /** Fitter for y. */
042 private final PolynomialFitter yFitter;
043
044 /**
045 * Default constructor. The degree of the fitting polynomials is set to 3.
046 */
047 public SmoothingPolynomialBicubicSplineInterpolator() {
048 this(3);
049 }
050
051 /**
052 * @param degree Degree of the polynomial fitting functions.
053 */
054 public SmoothingPolynomialBicubicSplineInterpolator(int degree) {
055 this(degree, degree);
056 }
057
058 /**
059 * @param xDegree Degree of the polynomial fitting functions along the
060 * x-dimension.
061 * @param yDegree Degree of the polynomial fitting functions along the
062 * y-dimension.
063 */
064 public SmoothingPolynomialBicubicSplineInterpolator(int xDegree,
065 int yDegree) {
066 xFitter = new PolynomialFitter(xDegree, new GaussNewtonOptimizer(false));
067 yFitter = new PolynomialFitter(yDegree, new GaussNewtonOptimizer(false));
068 }
069
070 /**
071 * {@inheritDoc}
072 */
073 @Override
074 public BicubicSplineInterpolatingFunction interpolate(final double[] xval,
075 final double[] yval,
076 final double[][] fval)
077 throws MathException {
078 if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
079 throw new NoDataException();
080 }
081 if (xval.length != fval.length) {
082 throw new DimensionMismatchException(xval.length, fval.length);
083 }
084
085 final int xLen = xval.length;
086 final int yLen = yval.length;
087
088 for (int i = 0; i < xLen; i++) {
089 if (fval[i].length != yLen) {
090 throw new DimensionMismatchException(fval[i].length, yLen);
091 }
092 }
093
094 MathUtils.checkOrder(xval);
095 MathUtils.checkOrder(yval);
096
097 // For each line y[j] (0 <= j < yLen), construct a polynomial, with
098 // respect to variable x, fitting array fval[][j]
099 final PolynomialFunction[] yPolyX = new PolynomialFunction[yLen];
100 for (int j = 0; j < yLen; j++) {
101 xFitter.clearObservations();
102 for (int i = 0; i < xLen; i++) {
103 xFitter.addObservedPoint(1, xval[i], fval[i][j]);
104 }
105
106 yPolyX[j] = xFitter.fit();
107 }
108
109 // For every knot (xval[i], yval[j]) of the grid, calculate corrected
110 // values fval_1
111 final double[][] fval_1 = new double[xLen][yLen];
112 for (int j = 0; j < yLen; j++) {
113 final PolynomialFunction f = yPolyX[j];
114 for (int i = 0; i < xLen; i++) {
115 fval_1[i][j] = f.value(xval[i]);
116 }
117 }
118
119 // For each line x[i] (0 <= i < xLen), construct a polynomial, with
120 // respect to variable y, fitting array fval_1[i][]
121 final PolynomialFunction[] xPolyY = new PolynomialFunction[xLen];
122 for (int i = 0; i < xLen; i++) {
123 yFitter.clearObservations();
124 for (int j = 0; j < yLen; j++) {
125 yFitter.addObservedPoint(1, yval[j], fval_1[i][j]);
126 }
127
128 xPolyY[i] = yFitter.fit();
129 }
130
131 // For every knot (xval[i], yval[j]) of the grid, calculate corrected
132 // values fval_2
133 final double[][] fval_2 = new double[xLen][yLen];
134 for (int i = 0; i < xLen; i++) {
135 final PolynomialFunction f = xPolyY[i];
136 for (int j = 0; j < yLen; j++) {
137 fval_2[i][j] = f.value(yval[j]);
138 }
139 }
140
141 return super.interpolate(xval, yval, fval_2);
142 }
143 }