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.DimensionMismatchException;
20import org.apache.commons.math.MathException;
21import org.apache.commons.math.analysis.UnivariateRealFunction;
22import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
23import org.apache.commons.math.exception.NoDataException;
24import org.apache.commons.math.util.MathUtils;
25
26/**
27 * Generates a bicubic interpolating function.
28 *
29 * @version $Revision: 980944 $ $Date: 2010-07-30 22:31:11 +0200 (ven. 30 juil. 2010) $
30 * @since 2.2
31 */
32public class BicubicSplineInterpolator
33    implements BivariateRealGridInterpolator {
34    /**
35     * {@inheritDoc}
36     */
37    public BicubicSplineInterpolatingFunction interpolate(final double[] xval,
38                                                          final double[] yval,
39                                                          final double[][] fval)
40        throws MathException, IllegalArgumentException {
41        if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
42            throw new NoDataException();
43        }
44        if (xval.length != fval.length) {
45            throw new DimensionMismatchException(xval.length, fval.length);
46        }
47
48        MathUtils.checkOrder(xval);
49        MathUtils.checkOrder(yval);
50
51        final int xLen = xval.length;
52        final int yLen = yval.length;
53
54        // Samples (first index is y-coordinate, i.e. subarray variable is x)
55        // 0 <= i < xval.length
56        // 0 <= j < yval.length
57        // fX[j][i] = f(xval[i], yval[j])
58        final double[][] fX = new double[yLen][xLen];
59        for (int i = 0; i < xLen; i++) {
60            if (fval[i].length != yLen) {
61                throw new DimensionMismatchException(fval[i].length, yLen);
62            }
63
64            for (int j = 0; j < yLen; j++) {
65                fX[j][i] = fval[i][j];
66            }
67        }
68
69        final SplineInterpolator spInterpolator = new SplineInterpolator();
70
71        // For each line y[j] (0 <= j < yLen), construct a 1D spline with
72        // respect to variable x
73        final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen];
74        for (int j = 0; j < yLen; j++) {
75            ySplineX[j] = spInterpolator.interpolate(xval, fX[j]);
76        }
77
78        // For each line x[i] (0 <= i < xLen), construct a 1D spline with
79        // respect to variable y generated by array fY_1[i]
80        final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen];
81        for (int i = 0; i < xLen; i++) {
82            xSplineY[i] = spInterpolator.interpolate(yval, fval[i]);
83        }
84
85        // Partial derivatives with respect to x at the grid knots
86        final double[][] dFdX = new double[xLen][yLen];
87        for (int j = 0; j < yLen; j++) {
88            final UnivariateRealFunction f = ySplineX[j].derivative();
89            for (int i = 0; i < xLen; i++) {
90                dFdX[i][j] = f.value(xval[i]);
91            }
92        }
93
94        // Partial derivatives with respect to y at the grid knots
95        final double[][] dFdY = new double[xLen][yLen];
96        for (int i = 0; i < xLen; i++) {
97            final UnivariateRealFunction f = xSplineY[i].derivative();
98            for (int j = 0; j < yLen; j++) {
99                dFdY[i][j] = f.value(yval[j]);
100            }
101        }
102
103        // Cross partial derivatives
104        final double[][] d2FdXdY = new double[xLen][yLen];
105        for (int i = 0; i < xLen ; i++) {
106            final int nI = nextIndex(i, xLen);
107            final int pI = previousIndex(i);
108            for (int j = 0; j < yLen; j++) {
109                final int nJ = nextIndex(j, yLen);
110                final int pJ = previousIndex(j);
111                d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] -
112                                 fval[pI][nJ] + fval[pI][pJ]) /
113                    ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]));
114            }
115        }
116
117        // Create the interpolating splines
118        return new BicubicSplineInterpolatingFunction(xval, yval, fval,
119                                                      dFdX, dFdY, d2FdXdY);
120    }
121
122    /**
123     * Compute the next index of an array, clipping if necessary.
124     * It is assumed (but not checked) that {@code i} is larger than or equal to 0}.
125     *
126     * @param i Index
127     * @param max Upper limit of the array
128     * @return the next index
129     */
130    private int nextIndex(int i, int max) {
131        final int index = i + 1;
132        return index < max ? index : index - 1;
133    }
134    /**
135     * Compute the previous index of an array, clipping if necessary.
136     * It is assumed (but not checked) that {@code i} is smaller than the size of the array.
137     *
138     * @param i Index
139     * @return the previous index
140     */
141    private int previousIndex(int i) {
142        final int index = i - 1;
143        return index >= 0 ? index : 0;
144    }
145}
146