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.math.analysis.interpolation; 018 019 import org.apache.commons.math.DimensionMismatchException; 020 import org.apache.commons.math.MathRuntimeException; 021 import org.apache.commons.math.MathException; 022 import org.apache.commons.math.util.MathUtils; 023 import org.apache.commons.math.analysis.UnivariateRealFunction; 024 import org.apache.commons.math.analysis.BivariateRealFunction; 025 import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; 026 027 /** 028 * Generates a bicubic interpolation function. 029 * Before interpolating, smoothing of the input data is performed using 030 * splines. 031 * See <b>Handbook on splines for the user</b>, ISBN 084939404X, 032 * chapter 2. 033 * 034 * @version $Revision$ $Date$ 035 * @since 2.1 036 */ 037 public class SmoothingBicubicSplineInterpolator 038 implements BivariateRealGridInterpolator { 039 /** 040 * {@inheritDoc} 041 */ 042 public BivariateRealFunction interpolate(final double[] xval, 043 final double[] yval, 044 final double[][] zval) 045 throws MathException, IllegalArgumentException { 046 if (xval.length == 0 || yval.length == 0 || zval.length == 0) { 047 throw MathRuntimeException.createIllegalArgumentException("no data"); 048 } 049 if (xval.length != zval.length) { 050 throw new DimensionMismatchException(xval.length, zval.length); 051 } 052 053 MathUtils.checkOrder(xval, 1, true); 054 MathUtils.checkOrder(yval, 1, true); 055 056 final int xLen = xval.length; 057 final int yLen = yval.length; 058 059 // Samples (first index is y-coordinate, i.e. subarray variable is x) 060 // 0 <= i < xval.length 061 // 0 <= j < yval.length 062 // zX[j][i] = f(xval[i], yval[j]) 063 final double[][] zX = new double[yLen][xLen]; 064 for (int i = 0; i < xLen; i++) { 065 if (zval[i].length != yLen) { 066 throw new DimensionMismatchException(zval[i].length, yLen); 067 } 068 069 for (int j = 0; j < yLen; j++) { 070 zX[j][i] = zval[i][j]; 071 } 072 } 073 074 final SplineInterpolator spInterpolator = new SplineInterpolator(); 075 076 // For each line y[j] (0 <= j < yLen), construct a 1D spline with 077 // respect to variable x 078 final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen]; 079 for (int j = 0; j < yLen; j++) { 080 ySplineX[j] = spInterpolator.interpolate(xval, zX[j]); 081 } 082 083 // For every knot (xval[i], yval[j]) of the grid, calculate corrected 084 // values zY_1 085 final double[][] zY_1 = new double[xLen][yLen]; 086 for (int j = 0; j < yLen; j++) { 087 final PolynomialSplineFunction f = ySplineX[j]; 088 for (int i = 0; i < xLen; i++) { 089 zY_1[i][j] = f.value(xval[i]); 090 } 091 } 092 093 // For each line x[i] (0 <= i < xLen), construct a 1D spline with 094 // respect to variable y generated by array zY_1[i] 095 final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen]; 096 for (int i = 0; i < xLen; i++) { 097 xSplineY[i] = spInterpolator.interpolate(yval, zY_1[i]); 098 } 099 100 // For every knot (xval[i], yval[j]) of the grid, calculate corrected 101 // values zY_2 102 final double[][] zY_2 = new double[xLen][yLen]; 103 for (int i = 0; i < xLen; i++) { 104 final PolynomialSplineFunction f = xSplineY[i]; 105 for (int j = 0; j < yLen; j++) { 106 zY_2[i][j] = f.value(yval[j]); 107 } 108 } 109 110 // Partial derivatives with respect to x at the grid knots 111 final double[][] dZdX = new double[xLen][yLen]; 112 for (int j = 0; j < yLen; j++) { 113 final UnivariateRealFunction f = ySplineX[j].derivative(); 114 for (int i = 0; i < xLen; i++) { 115 dZdX[i][j] = f.value(xval[i]); 116 } 117 } 118 119 // Partial derivatives with respect to y at the grid knots 120 final double[][] dZdY = new double[xLen][yLen]; 121 for (int i = 0; i < xLen; i++) { 122 final UnivariateRealFunction f = xSplineY[i].derivative(); 123 for (int j = 0; j < yLen; j++) { 124 dZdY[i][j] = f.value(yval[j]); 125 } 126 } 127 128 // Cross partial derivatives 129 final double[][] dZdXdY = new double[xLen][yLen]; 130 for (int i = 0; i < xLen ; i++) { 131 final int nI = nextIndex(i, xLen); 132 final int pI = previousIndex(i); 133 for (int j = 0; j < yLen; j++) { 134 final int nJ = nextIndex(j, yLen); 135 final int pJ = previousIndex(j); 136 dZdXdY[i][j] = (zY_2[nI][nJ] - zY_2[nI][pJ] - 137 zY_2[pI][nJ] + zY_2[pI][pJ]) / 138 ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ])) ; 139 } 140 } 141 142 // Create the interpolating splines 143 return new BicubicSplineInterpolatingFunction(xval, yval, zY_2, 144 dZdX, dZdY, dZdXdY); 145 } 146 147 /** 148 * Compute the next index of an array, clipping if necessary. 149 * It is assumed (but not checked) that {@code i} is larger than or equal to 0}. 150 * 151 * @param i Index 152 * @param max Upper limit of the array 153 * @return the next index 154 */ 155 private int nextIndex(int i, int max) { 156 final int index = i + 1; 157 return index < max ? index : index - 1; 158 } 159 /** 160 * Compute the previous index of an array, clipping if necessary. 161 * It is assumed (but not checked) that {@code i} is smaller than the size of the array. 162 * 163 * @param i Index 164 * @return the previous index 165 */ 166 private int previousIndex(int i) { 167 final int index = i - 1; 168 return index >= 0 ? index : 0; 169 } 170 }