Home | History | Annotate | Download | only in interpolation 
      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  */
     17 package org.apache.commons.math.analysis.interpolation;
     18 
     19 import org.apache.commons.math.exception.DimensionMismatchException;
     20 import org.apache.commons.math.exception.NoDataException;
     21 import org.apache.commons.math.MathException;
     22 import org.apache.commons.math.util.MathUtils;
     23 import org.apache.commons.math.optimization.general.GaussNewtonOptimizer;
     24 import org.apache.commons.math.optimization.fitting.PolynomialFitter;
     25 import 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  */
     35 public 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