1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17package org.apache.commons.math.analysis.interpolation;
18
19import org.apache.commons.math.exception.DimensionMismatchException;
20import org.apache.commons.math.exception.NoDataException;
21import org.apache.commons.math.MathException;
22import org.apache.commons.math.util.MathUtils;
23import org.apache.commons.math.optimization.general.GaussNewtonOptimizer;
24import org.apache.commons.math.optimization.fitting.PolynomialFitter;
25import org.apache.commons.math.analysis.polynomials.PolynomialFunction;
26
27/**
28 * Generates a bicubic interpolation function.
29 * Prior to generating the interpolating function, the input is smoothed using
30 * polynomial fitting.
31 *
32 * @version $Revision: 1003892 $ $Date: 2010-10-02 23:28:56 +0200 (sam. 02 oct. 2010) $
33 * @since 2.2
34 */
35public class SmoothingPolynomialBicubicSplineInterpolator
36    extends BicubicSplineInterpolator {
37
38    /** Fitter for x. */
39    private final PolynomialFitter xFitter;
40
41    /** Fitter for y. */
42    private final PolynomialFitter yFitter;
43
44    /**
45     * Default constructor. The degree of the fitting polynomials is set to 3.
46     */
47    public SmoothingPolynomialBicubicSplineInterpolator() {
48        this(3);
49    }
50
51    /**
52     * @param degree Degree of the polynomial fitting functions.
53     */
54    public SmoothingPolynomialBicubicSplineInterpolator(int degree) {
55        this(degree, degree);
56    }
57
58    /**
59     * @param xDegree Degree of the polynomial fitting functions along the
60     * x-dimension.
61     * @param yDegree Degree of the polynomial fitting functions along the
62     * y-dimension.
63     */
64    public SmoothingPolynomialBicubicSplineInterpolator(int xDegree,
65                                                        int yDegree) {
66        xFitter = new PolynomialFitter(xDegree, new GaussNewtonOptimizer(false));
67        yFitter = new PolynomialFitter(yDegree, new GaussNewtonOptimizer(false));
68    }
69
70    /**
71     * {@inheritDoc}
72     */
73    @Override
74    public BicubicSplineInterpolatingFunction interpolate(final double[] xval,
75                                                          final double[] yval,
76                                                          final double[][] fval)
77        throws MathException {
78        if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
79            throw new NoDataException();
80        }
81        if (xval.length != fval.length) {
82            throw new DimensionMismatchException(xval.length, fval.length);
83        }
84
85        final int xLen = xval.length;
86        final int yLen = yval.length;
87
88        for (int i = 0; i < xLen; i++) {
89            if (fval[i].length != yLen) {
90                throw new DimensionMismatchException(fval[i].length, yLen);
91            }
92        }
93
94        MathUtils.checkOrder(xval);
95        MathUtils.checkOrder(yval);
96
97        // For each line y[j] (0 <= j < yLen), construct a polynomial, with
98        // respect to variable x, fitting array fval[][j]
99        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}
144