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