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.math3.analysis.interpolation;
018
019 import org.apache.commons.math3.exception.DimensionMismatchException;
020 import org.apache.commons.math3.exception.NoDataException;
021 import org.apache.commons.math3.exception.NotPositiveException;
022 import org.apache.commons.math3.util.MathArrays;
023 import org.apache.commons.math3.util.Precision;
024 import org.apache.commons.math3.optim.nonlinear.vector.jacobian.GaussNewtonOptimizer;
025 import org.apache.commons.math3.fitting.PolynomialFitter;
026 import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
027 import org.apache.commons.math3.optim.SimpleVectorValueChecker;
028
029 /**
030 * Generates a bicubic interpolation function.
031 * Prior to generating the interpolating function, the input is smoothed using
032 * polynomial fitting.
033 *
034 * @version $Id: SmoothingPolynomialBicubicSplineInterpolator.java 1423558 2012-12-18 18:07:45Z erans $
035 * @since 2.2
036 */
037 public class SmoothingPolynomialBicubicSplineInterpolator
038 extends BicubicSplineInterpolator {
039 /** Fitter for x. */
040 private final PolynomialFitter xFitter;
041 /** Degree of the fitting polynomial. */
042 private final int xDegree;
043 /** Fitter for y. */
044 private final PolynomialFitter yFitter;
045 /** Degree of the fitting polynomial. */
046 private final int yDegree;
047
048 /**
049 * Default constructor. The degree of the fitting polynomials is set to 3.
050 */
051 public SmoothingPolynomialBicubicSplineInterpolator() {
052 this(3);
053 }
054
055 /**
056 * @param degree Degree of the polynomial fitting functions.
057 */
058 public SmoothingPolynomialBicubicSplineInterpolator(int degree) {
059 this(degree, degree);
060 }
061
062 /**
063 * @param xDegree Degree of the polynomial fitting functions along the
064 * x-dimension.
065 * @param yDegree Degree of the polynomial fitting functions along the
066 * y-dimension.
067 */
068 public SmoothingPolynomialBicubicSplineInterpolator(int xDegree,
069 int yDegree) {
070 if (xDegree < 0) {
071 throw new NotPositiveException(xDegree);
072 }
073 if (yDegree < 0) {
074 throw new NotPositiveException(yDegree);
075 }
076 this.xDegree = xDegree;
077 this.yDegree = yDegree;
078
079 final double safeFactor = 1e2;
080 final SimpleVectorValueChecker checker
081 = new SimpleVectorValueChecker(safeFactor * Precision.EPSILON,
082 safeFactor * Precision.SAFE_MIN);
083 xFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker));
084 yFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker));
085 }
086
087 /**
088 * {@inheritDoc}
089 */
090 @Override
091 public BicubicSplineInterpolatingFunction interpolate(final double[] xval,
092 final double[] yval,
093 final double[][] fval)
094 throws NoDataException,
095 DimensionMismatchException {
096 if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
097 throw new NoDataException();
098 }
099 if (xval.length != fval.length) {
100 throw new DimensionMismatchException(xval.length, fval.length);
101 }
102
103 final int xLen = xval.length;
104 final int yLen = yval.length;
105
106 for (int i = 0; i < xLen; i++) {
107 if (fval[i].length != yLen) {
108 throw new DimensionMismatchException(fval[i].length, yLen);
109 }
110 }
111
112 MathArrays.checkOrder(xval);
113 MathArrays.checkOrder(yval);
114
115 // For each line y[j] (0 <= j < yLen), construct a polynomial, with
116 // respect to variable x, fitting array fval[][j]
117 final PolynomialFunction[] yPolyX = new PolynomialFunction[yLen];
118 for (int j = 0; j < yLen; j++) {
119 xFitter.clearObservations();
120 for (int i = 0; i < xLen; i++) {
121 xFitter.addObservedPoint(1, xval[i], fval[i][j]);
122 }
123
124 // Initial guess for the fit is zero for each coefficients (of which
125 // there are "xDegree" + 1).
126 yPolyX[j] = new PolynomialFunction(xFitter.fit(new double[xDegree + 1]));
127 }
128
129 // For every knot (xval[i], yval[j]) of the grid, calculate corrected
130 // values fval_1
131 final double[][] fval_1 = new double[xLen][yLen];
132 for (int j = 0; j < yLen; j++) {
133 final PolynomialFunction f = yPolyX[j];
134 for (int i = 0; i < xLen; i++) {
135 fval_1[i][j] = f.value(xval[i]);
136 }
137 }
138
139 // For each line x[i] (0 <= i < xLen), construct a polynomial, with
140 // respect to variable y, fitting array fval_1[i][]
141 final PolynomialFunction[] xPolyY = new PolynomialFunction[xLen];
142 for (int i = 0; i < xLen; i++) {
143 yFitter.clearObservations();
144 for (int j = 0; j < yLen; j++) {
145 yFitter.addObservedPoint(1, yval[j], fval_1[i][j]);
146 }
147
148 // Initial guess for the fit is zero for each coefficients (of which
149 // there are "yDegree" + 1).
150 xPolyY[i] = new PolynomialFunction(yFitter.fit(new double[yDegree + 1]));
151 }
152
153 // For every knot (xval[i], yval[j]) of the grid, calculate corrected
154 // values fval_2
155 final double[][] fval_2 = new double[xLen][yLen];
156 for (int i = 0; i < xLen; i++) {
157 final PolynomialFunction f = xPolyY[i];
158 for (int j = 0; j < yLen; j++) {
159 fval_2[i][j] = f.value(yval[j]);
160 }
161 }
162
163 return super.interpolate(xval, yval, fval_2);
164 }
165 }