http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolatingFunction.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolatingFunction.java deleted file mode 100644 index b747841..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolatingFunction.java +++ /dev/null @@ -1,250 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.interpolation; - -import java.util.ArrayList; -import java.util.HashMap; -import java.util.List; -import java.util.Map; - -import org.apache.commons.math3.analysis.MultivariateFunction; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.linear.ArrayRealVector; -import org.apache.commons.math3.linear.RealVector; -import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator; -import org.apache.commons.math3.util.FastMath; - -/** - * Interpolating function that implements the - * <a href="http://www.dudziak.com/microsphere.php";>Microsphere Projection</a>. - * - */ -public class MicrosphereInterpolatingFunction - implements MultivariateFunction { - /** - * Space dimension. - */ - private final int dimension; - /** - * Internal accounting data for the interpolation algorithm. - * Each element of the list corresponds to one surface element of - * the microsphere. - */ - private final List<MicrosphereSurfaceElement> microsphere; - /** - * Exponent used in the power law that computes the weights of the - * sample data. - */ - private final double brightnessExponent; - /** - * Sample data. - */ - private final Map<RealVector, Double> samples; - - /** - * Class for storing the accounting data needed to perform the - * microsphere projection. - */ - private static class MicrosphereSurfaceElement { - /** Normal vector characterizing a surface element. */ - private final RealVector normal; - /** Illumination received from the brightest sample. */ - private double brightestIllumination; - /** Brightest sample. */ - private Map.Entry<RealVector, Double> brightestSample; - - /** - * @param n Normal vector characterizing a surface element - * of the microsphere. - */ - MicrosphereSurfaceElement(double[] n) { - normal = new ArrayRealVector(n); - } - - /** - * Return the normal vector. - * @return the normal vector - */ - RealVector normal() { - return normal; - } - - /** - * Reset "illumination" and "sampleIndex". - */ - void reset() { - brightestIllumination = 0; - brightestSample = null; - } - - /** - * Store the illumination and index of the brightest sample. - * @param illuminationFromSample illumination received from sample - * @param sample current sample illuminating the element - */ - void store(final double illuminationFromSample, - final Map.Entry<RealVector, Double> sample) { - if (illuminationFromSample > this.brightestIllumination) { - this.brightestIllumination = illuminationFromSample; - this.brightestSample = sample; - } - } - - /** - * Get the illumination of the element. - * @return the illumination. - */ - double illumination() { - return brightestIllumination; - } - - /** - * Get the sample illuminating the element the most. - * @return the sample. - */ - Map.Entry<RealVector, Double> sample() { - return brightestSample; - } - } - - /** - * @param xval Arguments for the interpolation points. - * {@code xval[i][0]} is the first component of interpolation point - * {@code i}, {@code xval[i][1]} is the second component, and so on - * until {@code xval[i][d-1]}, the last component of that interpolation - * point (where {@code dimension} is thus the dimension of the sampled - * space). - * @param yval Values for the interpolation points. - * @param brightnessExponent Brightness dimming factor. - * @param microsphereElements Number of surface elements of the - * microsphere. - * @param rand Unit vector generator for creating the microsphere. - * @throws DimensionMismatchException if the lengths of {@code yval} and - * {@code xval} (equal to {@code n}, the number of interpolation points) - * do not match, or the the arrays {@code xval[0]} ... {@code xval[n]}, - * have lengths different from {@code dimension}. - * @throws NoDataException if there an array has zero-length. - * @throws NullArgumentException if an argument is {@code null}. - */ - public MicrosphereInterpolatingFunction(double[][] xval, - double[] yval, - int brightnessExponent, - int microsphereElements, - UnitSphereRandomVectorGenerator rand) - throws DimensionMismatchException, - NoDataException, - NullArgumentException { - if (xval == null || - yval == null) { - throw new NullArgumentException(); - } - if (xval.length == 0) { - throw new NoDataException(); - } - if (xval.length != yval.length) { - throw new DimensionMismatchException(xval.length, yval.length); - } - if (xval[0] == null) { - throw new NullArgumentException(); - } - - dimension = xval[0].length; - this.brightnessExponent = brightnessExponent; - - // Copy data samples. - samples = new HashMap<RealVector, Double>(yval.length); - for (int i = 0; i < xval.length; ++i) { - final double[] xvalI = xval[i]; - if (xvalI == null) { - throw new NullArgumentException(); - } - if (xvalI.length != dimension) { - throw new DimensionMismatchException(xvalI.length, dimension); - } - - samples.put(new ArrayRealVector(xvalI), yval[i]); - } - - microsphere = new ArrayList<MicrosphereSurfaceElement>(microsphereElements); - // Generate the microsphere, assuming that a fairly large number of - // randomly generated normals will represent a sphere. - for (int i = 0; i < microsphereElements; i++) { - microsphere.add(new MicrosphereSurfaceElement(rand.nextVector())); - } - } - - /** - * @param point Interpolation point. - * @return the interpolated value. - * @throws DimensionMismatchException if point dimension does not math sample - */ - public double value(double[] point) throws DimensionMismatchException { - final RealVector p = new ArrayRealVector(point); - - // Reset. - for (MicrosphereSurfaceElement md : microsphere) { - md.reset(); - } - - // Compute contribution of each sample points to the microsphere elements illumination - for (Map.Entry<RealVector, Double> sd : samples.entrySet()) { - - // Vector between interpolation point and current sample point. - final RealVector diff = sd.getKey().subtract(p); - final double diffNorm = diff.getNorm(); - - if (FastMath.abs(diffNorm) < FastMath.ulp(1d)) { - // No need to interpolate, as the interpolation point is - // actually (very close to) one of the sampled points. - return sd.getValue(); - } - - for (MicrosphereSurfaceElement md : microsphere) { - final double w = FastMath.pow(diffNorm, -brightnessExponent); - md.store(cosAngle(diff, md.normal()) * w, sd); - } - - } - - // Interpolation calculation. - double value = 0; - double totalWeight = 0; - for (MicrosphereSurfaceElement md : microsphere) { - final double iV = md.illumination(); - final Map.Entry<RealVector, Double> sd = md.sample(); - if (sd != null) { - value += iV * sd.getValue(); - totalWeight += iV; - } - } - - return value / totalWeight; - } - - /** - * Compute the cosine of the angle between 2 vectors. - * - * @param v Vector. - * @param w Vector. - * @return the cosine of the angle between {@code v} and {@code w}. - */ - private double cosAngle(final RealVector v, final RealVector w) { - return v.dotProduct(w) / (v.getNorm() * w.getNorm()); - } -}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolator.java deleted file mode 100644 index c9881ce..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/MicrosphereInterpolator.java +++ /dev/null @@ -1,102 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.interpolation; - -import org.apache.commons.math3.analysis.MultivariateFunction; -import org.apache.commons.math3.exception.NotPositiveException; -import org.apache.commons.math3.exception.NotStrictlyPositiveException; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.random.UnitSphereRandomVectorGenerator; - -/** - * Interpolator that implements the algorithm described in - * <em>William Dudziak</em>'s - * <a href="http://www.dudziak.com/microsphere.pdf";>MS thesis</a>. - * - * @since 2.1 - */ -public class MicrosphereInterpolator - implements MultivariateInterpolator { - /** - * Default number of surface elements that composes the microsphere. - */ - public static final int DEFAULT_MICROSPHERE_ELEMENTS = 2000; - /** - * Default exponent used the weights calculation. - */ - public static final int DEFAULT_BRIGHTNESS_EXPONENT = 2; - /** - * Number of surface elements of the microsphere. - */ - private final int microsphereElements; - /** - * Exponent used in the power law that computes the weights of the - * sample data. - */ - private final int brightnessExponent; - - /** - * Create a microsphere interpolator with default settings. - * Calling this constructor is equivalent to call {@link - * #MicrosphereInterpolator(int, int) - * MicrosphereInterpolator(MicrosphereInterpolator.DEFAULT_MICROSPHERE_ELEMENTS, - * MicrosphereInterpolator.DEFAULT_BRIGHTNESS_EXPONENT)}. - */ - public MicrosphereInterpolator() { - this(DEFAULT_MICROSPHERE_ELEMENTS, DEFAULT_BRIGHTNESS_EXPONENT); - } - - /** Create a microsphere interpolator. - * @param elements Number of surface elements of the microsphere. - * @param exponent Exponent used in the power law that computes the - * weights (distance dimming factor) of the sample data. - * @throws NotPositiveException if {@code exponent < 0}. - * @throws NotStrictlyPositiveException if {@code elements <= 0}. - */ - public MicrosphereInterpolator(final int elements, - final int exponent) - throws NotPositiveException, - NotStrictlyPositiveException { - if (exponent < 0) { - throw new NotPositiveException(exponent); - } - if (elements <= 0) { - throw new NotStrictlyPositiveException(elements); - } - - microsphereElements = elements; - brightnessExponent = exponent; - } - - /** - * {@inheritDoc} - */ - public MultivariateFunction interpolate(final double[][] xval, - final double[] yval) - throws DimensionMismatchException, - NoDataException, - NullArgumentException { - final UnitSphereRandomVectorGenerator rand - = new UnitSphereRandomVectorGenerator(xval[0].length); - return new MicrosphereInterpolatingFunction(xval, yval, - brightnessExponent, - microsphereElements, - rand); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/MultivariateInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/MultivariateInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/MultivariateInterpolator.java deleted file mode 100644 index 7d76374..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/MultivariateInterpolator.java +++ /dev/null @@ -1,51 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.interpolation; - -import org.apache.commons.math3.analysis.MultivariateFunction; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.MathIllegalArgumentException; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.NullArgumentException; - -/** - * Interface representing a univariate real interpolating function. - * - * @since 2.1 - */ -public interface MultivariateInterpolator { - - /** - * Computes an interpolating function for the data set. - * - * @param xval the arguments for the interpolation points. - * {@code xval[i][0]} is the first component of interpolation point - * {@code i}, {@code xval[i][1]} is the second component, and so on - * until {@code xval[i][d-1]}, the last component of that interpolation - * point (where {@code d} is thus the dimension of the space). - * @param yval the values for the interpolation points - * @return a function which interpolates the data set - * @throws MathIllegalArgumentException if the arguments violate assumptions - * made by the interpolation algorithm. - * @throws DimensionMismatchException when the array dimensions are not consistent. - * @throws NoDataException if an array has zero-length. - * @throws NullArgumentException if the arguments are {@code null}. - */ - MultivariateFunction interpolate(double[][] xval, double[] yval) - throws MathIllegalArgumentException, DimensionMismatchException, - NoDataException, NullArgumentException; -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/NevilleInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/NevilleInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/NevilleInterpolator.java deleted file mode 100644 index 6b47451..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/NevilleInterpolator.java +++ /dev/null @@ -1,60 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.interpolation; - -import java.io.Serializable; - -import org.apache.commons.math3.analysis.polynomials.PolynomialFunctionLagrangeForm; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NumberIsTooSmallException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; - -/** - * Implements the <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html";> - * Neville's Algorithm</a> for interpolation of real univariate functions. For - * reference, see <b>Introduction to Numerical Analysis</b>, ISBN 038795452X, - * chapter 2. - * <p> - * The actual code of Neville's algorithm is in PolynomialFunctionLagrangeForm, - * this class provides an easy-to-use interface to it.</p> - * - * @since 1.2 - */ -public class NevilleInterpolator implements UnivariateInterpolator, - Serializable { - - /** serializable version identifier */ - static final long serialVersionUID = 3003707660147873733L; - - /** - * Computes an interpolating function for the data set. - * - * @param x Interpolating points. - * @param y Interpolating values. - * @return a function which interpolates the data set - * @throws DimensionMismatchException if the array lengths are different. - * @throws NumberIsTooSmallException if the number of points is less than 2. - * @throws NonMonotonicSequenceException if two abscissae have the same - * value. - */ - public PolynomialFunctionLagrangeForm interpolate(double x[], double y[]) - throws DimensionMismatchException, - NumberIsTooSmallException, - NonMonotonicSequenceException { - return new PolynomialFunctionLagrangeForm(x, y); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunction.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunction.java deleted file mode 100644 index 7dd135a..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunction.java +++ /dev/null @@ -1,210 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.interpolation; - -import java.util.Arrays; -import org.apache.commons.math3.analysis.BivariateFunction; -import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.InsufficientDataException; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.exception.OutOfRangeException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.util.MathArrays; - -/** - * Function that implements the - * <a href="http://www.paulinternet.nl/?page=bicubic";>bicubic spline</a> - * interpolation. - * This implementation currently uses {@link AkimaSplineInterpolator} as the - * underlying one-dimensional interpolator, which requires 5 sample points; - * insufficient data will raise an exception when the - * {@link #value(double,double) value} method is called. - * - * @since 3.4 - */ -public class PiecewiseBicubicSplineInterpolatingFunction - implements BivariateFunction { - /** The minimum number of points that are needed to compute the function. */ - private static final int MIN_NUM_POINTS = 5; - /** Samples x-coordinates */ - private final double[] xval; - /** Samples y-coordinates */ - private final double[] yval; - /** Set of cubic splines patching the whole data grid */ - private final double[][] fval; - - /** - * @param x Sample values of the x-coordinate, in increasing order. - * @param y Sample values of the y-coordinate, in increasing order. - * @param f Values of the function on every grid point. the expected number - * of elements. - * @throws NonMonotonicSequenceException if {@code x} or {@code y} are not - * strictly increasing. - * @throws NullArgumentException if any of the arguments are null - * @throws NoDataException if any of the arrays has zero length. - * @throws DimensionMismatchException if the length of x and y don't match the row, column - * height of f - */ - public PiecewiseBicubicSplineInterpolatingFunction(double[] x, - double[] y, - double[][] f) - throws DimensionMismatchException, - NullArgumentException, - NoDataException, - NonMonotonicSequenceException { - if (x == null || - y == null || - f == null || - f[0] == null) { - throw new NullArgumentException(); - } - - final int xLen = x.length; - final int yLen = y.length; - - if (xLen == 0 || - yLen == 0 || - f.length == 0 || - f[0].length == 0) { - throw new NoDataException(); - } - - if (xLen < MIN_NUM_POINTS || - yLen < MIN_NUM_POINTS || - f.length < MIN_NUM_POINTS || - f[0].length < MIN_NUM_POINTS) { - throw new InsufficientDataException(); - } - - if (xLen != f.length) { - throw new DimensionMismatchException(xLen, f.length); - } - - if (yLen != f[0].length) { - throw new DimensionMismatchException(yLen, f[0].length); - } - - MathArrays.checkOrder(x); - MathArrays.checkOrder(y); - - xval = x.clone(); - yval = y.clone(); - fval = f.clone(); - } - - /** - * {@inheritDoc} - */ - public double value(double x, - double y) - throws OutOfRangeException { - final AkimaSplineInterpolator interpolator = new AkimaSplineInterpolator(); - final int offset = 2; - final int count = offset + 3; - final int i = searchIndex(x, xval, offset, count); - final int j = searchIndex(y, yval, offset, count); - - final double xArray[] = new double[count]; - final double yArray[] = new double[count]; - final double zArray[] = new double[count]; - final double interpArray[] = new double[count]; - - for (int index = 0; index < count; index++) { - xArray[index] = xval[i + index]; - yArray[index] = yval[j + index]; - } - - for (int zIndex = 0; zIndex < count; zIndex++) { - for (int index = 0; index < count; index++) { - zArray[index] = fval[i + index][j + zIndex]; - } - final PolynomialSplineFunction spline = interpolator.interpolate(xArray, zArray); - interpArray[zIndex] = spline.value(x); - } - - final PolynomialSplineFunction spline = interpolator.interpolate(yArray, interpArray); - - double returnValue = spline.value(y); - - return returnValue; - } - - /** - * Indicates whether a point is within the interpolation range. - * - * @param x First coordinate. - * @param y Second coordinate. - * @return {@code true} if (x, y) is a valid point. - * @since 3.3 - */ - public boolean isValidPoint(double x, - double y) { - if (x < xval[0] || - x > xval[xval.length - 1] || - y < yval[0] || - y > yval[yval.length - 1]) { - return false; - } else { - return true; - } - } - - /** - * @param c Coordinate. - * @param val Coordinate samples. - * @param offset how far back from found value to offset for querying - * @param count total number of elements forward from beginning that will be - * queried - * @return the index in {@code val} corresponding to the interval containing - * {@code c}. - * @throws OutOfRangeException if {@code c} is out of the range defined by - * the boundary values of {@code val}. - */ - private int searchIndex(double c, - double[] val, - int offset, - int count) { - int r = Arrays.binarySearch(val, c); - - if (r == -1 || r == -val.length - 1) { - throw new OutOfRangeException(c, val[0], val[val.length - 1]); - } - - if (r < 0) { - // "c" in within an interpolation sub-interval, which returns - // negative - // need to remove the negative sign for consistency - r = -r - offset - 1; - } else { - r -= offset; - } - - if (r < 0) { - r = 0; - } - - if ((r + count) >= val.length) { - // "c" is the last sample of the range: Return the index - // of the sample at the lower end of the last sub-interval. - r = val.length - count; - } - - return r; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolator.java deleted file mode 100644 index 826f328..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolator.java +++ /dev/null @@ -1,61 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.interpolation; - -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.util.MathArrays; - -/** - * Generates a piecewise-bicubic interpolating function. - * - * @since 2.2 - */ -public class PiecewiseBicubicSplineInterpolator - implements BivariateGridInterpolator { - - /** - * {@inheritDoc} - */ - public PiecewiseBicubicSplineInterpolatingFunction interpolate( final double[] xval, - final double[] yval, - final double[][] fval) - throws DimensionMismatchException, - NullArgumentException, - NoDataException, - NonMonotonicSequenceException { - if ( xval == null || - yval == null || - fval == null || - fval[0] == null ) { - throw new NullArgumentException(); - } - - if ( xval.length == 0 || - yval.length == 0 || - fval.length == 0 ) { - throw new NoDataException(); - } - - MathArrays.checkOrder(xval); - MathArrays.checkOrder(yval); - - return new PiecewiseBicubicSplineInterpolatingFunction( xval, yval, fval ); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java deleted file mode 100644 index e1639b2..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java +++ /dev/null @@ -1,171 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.interpolation; - -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.exception.NotPositiveException; -import org.apache.commons.math3.exception.NullArgumentException; -import org.apache.commons.math3.util.MathArrays; -import org.apache.commons.math3.util.Precision; -import org.apache.commons.math3.optim.nonlinear.vector.jacobian.GaussNewtonOptimizer; -import org.apache.commons.math3.fitting.PolynomialFitter; -import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; -import org.apache.commons.math3.optim.SimpleVectorValueChecker; - -/** - * Generates a bicubic interpolation function. - * Prior to generating the interpolating function, the input is smoothed using - * polynomial fitting. - * - * @since 2.2 - * @deprecated To be removed in 4.0 (see MATH-1166). - */ -@Deprecated -public class SmoothingPolynomialBicubicSplineInterpolator - extends BicubicSplineInterpolator { - /** Fitter for x. */ - private final PolynomialFitter xFitter; - /** Degree of the fitting polynomial. */ - private final int xDegree; - /** Fitter for y. */ - private final PolynomialFitter yFitter; - /** Degree of the fitting polynomial. */ - private final int yDegree; - - /** - * Default constructor. The degree of the fitting polynomials is set to 3. - */ - public SmoothingPolynomialBicubicSplineInterpolator() { - this(3); - } - - /** - * @param degree Degree of the polynomial fitting functions. - * @exception NotPositiveException if degree is not positive - */ - public SmoothingPolynomialBicubicSplineInterpolator(int degree) - throws NotPositiveException { - this(degree, degree); - } - - /** - * @param xDegree Degree of the polynomial fitting functions along the - * x-dimension. - * @param yDegree Degree of the polynomial fitting functions along the - * y-dimension. - * @exception NotPositiveException if degrees are not positive - */ - public SmoothingPolynomialBicubicSplineInterpolator(int xDegree, int yDegree) - throws NotPositiveException { - if (xDegree < 0) { - throw new NotPositiveException(xDegree); - } - if (yDegree < 0) { - throw new NotPositiveException(yDegree); - } - this.xDegree = xDegree; - this.yDegree = yDegree; - - final double safeFactor = 1e2; - final SimpleVectorValueChecker checker - = new SimpleVectorValueChecker(safeFactor * Precision.EPSILON, - safeFactor * Precision.SAFE_MIN); - xFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker)); - yFitter = new PolynomialFitter(new GaussNewtonOptimizer(false, checker)); - } - - /** - * {@inheritDoc} - */ - @Override - public BicubicSplineInterpolatingFunction interpolate(final double[] xval, - final double[] yval, - final double[][] fval) - throws NoDataException, NullArgumentException, - DimensionMismatchException, NonMonotonicSequenceException { - if (xval.length == 0 || yval.length == 0 || fval.length == 0) { - throw new NoDataException(); - } - if (xval.length != fval.length) { - throw new DimensionMismatchException(xval.length, fval.length); - } - - final int xLen = xval.length; - final int yLen = yval.length; - - for (int i = 0; i < xLen; i++) { - if (fval[i].length != yLen) { - throw new DimensionMismatchException(fval[i].length, yLen); - } - } - - MathArrays.checkOrder(xval); - MathArrays.checkOrder(yval); - - // For each line y[j] (0 <= j < yLen), construct a polynomial, with - // respect to variable x, fitting array fval[][j] - final PolynomialFunction[] yPolyX = new PolynomialFunction[yLen]; - for (int j = 0; j < yLen; j++) { - xFitter.clearObservations(); - for (int i = 0; i < xLen; i++) { - xFitter.addObservedPoint(1, xval[i], fval[i][j]); - } - - // Initial guess for the fit is zero for each coefficients (of which - // there are "xDegree" + 1). - yPolyX[j] = new PolynomialFunction(xFitter.fit(new double[xDegree + 1])); - } - - // For every knot (xval[i], yval[j]) of the grid, calculate corrected - // values fval_1 - final double[][] fval_1 = new double[xLen][yLen]; - for (int j = 0; j < yLen; j++) { - final PolynomialFunction f = yPolyX[j]; - for (int i = 0; i < xLen; i++) { - fval_1[i][j] = f.value(xval[i]); - } - } - - // For each line x[i] (0 <= i < xLen), construct a polynomial, with - // respect to variable y, fitting array fval_1[i][] - final PolynomialFunction[] xPolyY = new PolynomialFunction[xLen]; - for (int i = 0; i < xLen; i++) { - yFitter.clearObservations(); - for (int j = 0; j < yLen; j++) { - yFitter.addObservedPoint(1, yval[j], fval_1[i][j]); - } - - // Initial guess for the fit is zero for each coefficients (of which - // there are "yDegree" + 1). - xPolyY[i] = new PolynomialFunction(yFitter.fit(new double[yDegree + 1])); - } - - // For every knot (xval[i], yval[j]) of the grid, calculate corrected - // values fval_2 - final double[][] fval_2 = new double[xLen][yLen]; - for (int i = 0; i < xLen; i++) { - final PolynomialFunction f = xPolyY[i]; - for (int j = 0; j < yLen; j++) { - fval_2[i][j] = f.value(yval[j]); - } - } - - return super.interpolate(xval, yval, fval_2); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/SplineInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/SplineInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/SplineInterpolator.java deleted file mode 100644 index a9ca862..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/SplineInterpolator.java +++ /dev/null @@ -1,127 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.interpolation; - -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.util.LocalizedFormats; -import org.apache.commons.math3.exception.NumberIsTooSmallException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; -import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction; -import org.apache.commons.math3.util.MathArrays; - -/** - * Computes a natural (also known as "free", "unclamped") cubic spline interpolation for the data set. - * <p> - * The {@link #interpolate(double[], double[])} method returns a {@link PolynomialSplineFunction} - * consisting of n cubic polynomials, defined over the subintervals determined by the x values, - * x[0] < x[i] ... < x[n]. The x values are referred to as "knot points."</p> - * <p> - * The value of the PolynomialSplineFunction at a point x that is greater than or equal to the smallest - * knot point and strictly less than the largest knot point is computed by finding the subinterval to which - * x belongs and computing the value of the corresponding polynomial at <code>x - x[i] </code> where - * <code>i</code> is the index of the subinterval. See {@link PolynomialSplineFunction} for more details. - * </p> - * <p> - * The interpolating polynomials satisfy: <ol> - * <li>The value of the PolynomialSplineFunction at each of the input x values equals the - * corresponding y value.</li> - * <li>Adjacent polynomials are equal through two derivatives at the knot points (i.e., adjacent polynomials - * "match up" at the knot points, as do their first and second derivatives).</li> - * </ol></p> - * <p> - * The cubic spline interpolation algorithm implemented is as described in R.L. Burden, J.D. Faires, - * <u>Numerical Analysis</u>, 4th Ed., 1989, PWS-Kent, ISBN 0-53491-585-X, pp 126-131. - * </p> - * - */ -public class SplineInterpolator implements UnivariateInterpolator { - /** - * Computes an interpolating function for the data set. - * @param x the arguments for the interpolation points - * @param y the values for the interpolation points - * @return a function which interpolates the data set - * @throws DimensionMismatchException if {@code x} and {@code y} - * have different sizes. - * @throws NonMonotonicSequenceException if {@code x} is not sorted in - * strict increasing order. - * @throws NumberIsTooSmallException if the size of {@code x} is smaller - * than 3. - */ - public PolynomialSplineFunction interpolate(double x[], double y[]) - throws DimensionMismatchException, - NumberIsTooSmallException, - NonMonotonicSequenceException { - if (x.length != y.length) { - throw new DimensionMismatchException(x.length, y.length); - } - - if (x.length < 3) { - throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS, - x.length, 3, true); - } - - // Number of intervals. The number of data points is n + 1. - final int n = x.length - 1; - - MathArrays.checkOrder(x); - - // Differences between knot points - final double h[] = new double[n]; - for (int i = 0; i < n; i++) { - h[i] = x[i + 1] - x[i]; - } - - final double mu[] = new double[n]; - final double z[] = new double[n + 1]; - mu[0] = 0d; - z[0] = 0d; - double g = 0; - for (int i = 1; i < n; i++) { - g = 2d * (x[i+1] - x[i - 1]) - h[i - 1] * mu[i -1]; - mu[i] = h[i] / g; - z[i] = (3d * (y[i + 1] * h[i - 1] - y[i] * (x[i + 1] - x[i - 1])+ y[i - 1] * h[i]) / - (h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g; - } - - // cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants) - final double b[] = new double[n]; - final double c[] = new double[n + 1]; - final double d[] = new double[n]; - - z[n] = 0d; - c[n] = 0d; - - for (int j = n -1; j >=0; j--) { - c[j] = z[j] - mu[j] * c[j + 1]; - b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2d * c[j]) / 3d; - d[j] = (c[j + 1] - c[j]) / (3d * h[j]); - } - - final PolynomialFunction polynomials[] = new PolynomialFunction[n]; - final double coefficients[] = new double[4]; - for (int i = 0; i < n; i++) { - coefficients[0] = y[i]; - coefficients[1] = b[i]; - coefficients[2] = c[i]; - coefficients[3] = d[i]; - polynomials[i] = new PolynomialFunction(coefficients); - } - - return new PolynomialSplineFunction(x, polynomials); - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunction.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunction.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunction.java deleted file mode 100644 index 9344d89..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolatingFunction.java +++ /dev/null @@ -1,508 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.interpolation; - -import org.apache.commons.math3.analysis.TrivariateFunction; -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.OutOfRangeException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.util.MathArrays; - -/** - * Function that implements the - * <a href="http://en.wikipedia.org/wiki/Tricubic_interpolation";> - * tricubic spline interpolation</a>, as proposed in - * <blockquote> - * Tricubic interpolation in three dimensions<br> - * F. Lekien and J. Marsden<br> - * <em>Int. J. Numer. Meth. Eng</em> 2005; <b>63</b>:455-471<br> - * </blockquote> - * - * @since 3.4. - */ -public class TricubicInterpolatingFunction - implements TrivariateFunction { - /** - * Matrix to compute the spline coefficients from the function values - * and function derivatives values - */ - private static final double[][] AINV = { - { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { -3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { -3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 9,-9,-9,9,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,0,0,0,0,0,0,0,0,4,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { -6,6,6,-6,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,0,0,0,0,0,0,0,0,-2,-2,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { -6,6,6,-6,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 4,-4,-4,4,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,-9,9,0,0,0,0,0,0,0,0,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,4,2,2,1,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,-2,-2,-1,-1,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,-2,-1,-2,-1,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,-4,4,0,0,0,0,0,0,0,0,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,1,1,1,1,0,0,0,0 }, - {-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 9,-9,0,0,-9,9,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,0,0,0,0,0,0,0,0,4,2,0,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { -6,6,0,0,6,-6,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,0,0,0,0,0,0,0,0,-2,-2,0,0,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,0,0,-9,9,0,0,0,0,0,0,0,0,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,4,2,0,0,2,1,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,-2,-2,0,0,-1,-1,0,0 }, - { 9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0 }, - { -27,27,27,-27,27,-27,-27,27,-18,-9,18,9,18,9,-18,-9,-18,18,-9,9,18,-18,9,-9,-18,18,18,-18,-9,9,9,-9,-12,-6,-6,-3,12,6,6,3,-12,-6,12,6,-6,-3,6,3,-12,12,-6,6,-6,6,-3,3,-8,-4,-4,-2,-4,-2,-2,-1 }, - { 18,-18,-18,18,-18,18,18,-18,9,9,-9,-9,-9,-9,9,9,12,-12,6,-6,-12,12,-6,6,12,-12,-12,12,6,-6,-6,6,6,6,3,3,-6,-6,-3,-3,6,6,-6,-6,3,3,-3,-3,8,-8,4,-4,4,-4,2,-2,4,4,2,2,2,2,1,1 }, - { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0 }, - { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,9,-9,9,-9,-9,9,-9,9,12,-12,-12,12,6,-6,-6,6,6,3,6,3,-6,-3,-6,-3,8,4,-8,-4,4,2,-4,-2,6,-6,6,-6,3,-3,3,-3,4,2,4,2,2,1,2,1 }, - { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-6,6,-6,6,6,-6,6,-6,-8,8,8,-8,-4,4,4,-4,-3,-3,-3,-3,3,3,3,3,-4,-4,4,4,-2,-2,2,2,-4,4,-4,4,-2,2,-2,2,-2,-2,-2,-2,-1,-1,-1,-1 }, - { 2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { -6,6,0,0,6,-6,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 4,-4,0,0,-4,4,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,-2,-1,0,0,-2,-1,0,0 }, - { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,0,0,-4,4,0,0,0,0,0,0,0,0,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,1,1,0,0,1,1,0,0 }, - { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0 }, - { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,12,-12,6,-6,-12,12,-6,6,9,-9,-9,9,9,-9,-9,9,8,4,4,2,-8,-4,-4,-2,6,3,-6,-3,6,3,-6,-3,6,-6,3,-3,6,-6,3,-3,4,2,2,1,4,2,2,1 }, - { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-8,8,-4,4,8,-8,4,-4,-6,6,6,-6,-6,6,6,-6,-4,-4,-2,-2,4,4,2,2,-3,-3,3,3,-3,-3,3,3,-4,4,-2,2,-4,4,-2,2,-2,-2,-1,-1,-2,-2,-1,-1 }, - { 4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0 }, - { 0,0,0,0,0,0,0,0,4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0 }, - { -12,12,12,-12,12,-12,-12,12,-8,-4,8,4,8,4,-8,-4,-6,6,-6,6,6,-6,6,-6,-6,6,6,-6,-6,6,6,-6,-4,-2,-4,-2,4,2,4,2,-4,-2,4,2,-4,-2,4,2,-3,3,-3,3,-3,3,-3,3,-2,-1,-2,-1,-2,-1,-2,-1 }, - { 8,-8,-8,8,-8,8,8,-8,4,4,-4,-4,-4,-4,4,4,4,-4,4,-4,-4,4,-4,4,4,-4,-4,4,4,-4,-4,4,2,2,2,2,-2,-2,-2,-2,2,2,-2,-2,2,2,-2,-2,2,-2,2,-2,2,-2,2,-2,1,1,1,1,1,1,1,1 } - }; - - /** Samples x-coordinates */ - private final double[] xval; - /** Samples y-coordinates */ - private final double[] yval; - /** Samples z-coordinates */ - private final double[] zval; - /** Set of cubic splines pacthing the whole data grid */ - private final TricubicFunction[][][] splines; - - /** - * @param x Sample values of the x-coordinate, in increasing order. - * @param y Sample values of the y-coordinate, in increasing order. - * @param z Sample values of the y-coordinate, in increasing order. - * @param f Values of the function on every grid point. - * @param dFdX Values of the partial derivative of function with respect to x on every grid point. - * @param dFdY Values of the partial derivative of function with respect to y on every grid point. - * @param dFdZ Values of the partial derivative of function with respect to z on every grid point. - * @param d2FdXdY Values of the cross partial derivative of function on every grid point. - * @param d2FdXdZ Values of the cross partial derivative of function on every grid point. - * @param d2FdYdZ Values of the cross partial derivative of function on every grid point. - * @param d3FdXdYdZ Values of the cross partial derivative of function on every grid point. - * @throws NoDataException if any of the arrays has zero length. - * @throws DimensionMismatchException if the various arrays do not contain the expected number of elements. - * @throws NonMonotonicSequenceException if {@code x}, {@code y} or {@code z} are not strictly increasing. - */ - public TricubicInterpolatingFunction(double[] x, - double[] y, - double[] z, - double[][][] f, - double[][][] dFdX, - double[][][] dFdY, - double[][][] dFdZ, - double[][][] d2FdXdY, - double[][][] d2FdXdZ, - double[][][] d2FdYdZ, - double[][][] d3FdXdYdZ) - throws NoDataException, - DimensionMismatchException, - NonMonotonicSequenceException { - final int xLen = x.length; - final int yLen = y.length; - final int zLen = z.length; - - if (xLen == 0 || yLen == 0 || z.length == 0 || f.length == 0 || f[0].length == 0) { - throw new NoDataException(); - } - if (xLen != f.length) { - throw new DimensionMismatchException(xLen, f.length); - } - if (xLen != dFdX.length) { - throw new DimensionMismatchException(xLen, dFdX.length); - } - if (xLen != dFdY.length) { - throw new DimensionMismatchException(xLen, dFdY.length); - } - if (xLen != dFdZ.length) { - throw new DimensionMismatchException(xLen, dFdZ.length); - } - if (xLen != d2FdXdY.length) { - throw new DimensionMismatchException(xLen, d2FdXdY.length); - } - if (xLen != d2FdXdZ.length) { - throw new DimensionMismatchException(xLen, d2FdXdZ.length); - } - if (xLen != d2FdYdZ.length) { - throw new DimensionMismatchException(xLen, d2FdYdZ.length); - } - if (xLen != d3FdXdYdZ.length) { - throw new DimensionMismatchException(xLen, d3FdXdYdZ.length); - } - - MathArrays.checkOrder(x); - MathArrays.checkOrder(y); - MathArrays.checkOrder(z); - - xval = x.clone(); - yval = y.clone(); - zval = z.clone(); - - final int lastI = xLen - 1; - final int lastJ = yLen - 1; - final int lastK = zLen - 1; - splines = new TricubicFunction[lastI][lastJ][lastK]; - - for (int i = 0; i < lastI; i++) { - if (f[i].length != yLen) { - throw new DimensionMismatchException(f[i].length, yLen); - } - if (dFdX[i].length != yLen) { - throw new DimensionMismatchException(dFdX[i].length, yLen); - } - if (dFdY[i].length != yLen) { - throw new DimensionMismatchException(dFdY[i].length, yLen); - } - if (dFdZ[i].length != yLen) { - throw new DimensionMismatchException(dFdZ[i].length, yLen); - } - if (d2FdXdY[i].length != yLen) { - throw new DimensionMismatchException(d2FdXdY[i].length, yLen); - } - if (d2FdXdZ[i].length != yLen) { - throw new DimensionMismatchException(d2FdXdZ[i].length, yLen); - } - if (d2FdYdZ[i].length != yLen) { - throw new DimensionMismatchException(d2FdYdZ[i].length, yLen); - } - if (d3FdXdYdZ[i].length != yLen) { - throw new DimensionMismatchException(d3FdXdYdZ[i].length, yLen); - } - - final int ip1 = i + 1; - final double xR = xval[ip1] - xval[i]; - for (int j = 0; j < lastJ; j++) { - if (f[i][j].length != zLen) { - throw new DimensionMismatchException(f[i][j].length, zLen); - } - if (dFdX[i][j].length != zLen) { - throw new DimensionMismatchException(dFdX[i][j].length, zLen); - } - if (dFdY[i][j].length != zLen) { - throw new DimensionMismatchException(dFdY[i][j].length, zLen); - } - if (dFdZ[i][j].length != zLen) { - throw new DimensionMismatchException(dFdZ[i][j].length, zLen); - } - if (d2FdXdY[i][j].length != zLen) { - throw new DimensionMismatchException(d2FdXdY[i][j].length, zLen); - } - if (d2FdXdZ[i][j].length != zLen) { - throw new DimensionMismatchException(d2FdXdZ[i][j].length, zLen); - } - if (d2FdYdZ[i][j].length != zLen) { - throw new DimensionMismatchException(d2FdYdZ[i][j].length, zLen); - } - if (d3FdXdYdZ[i][j].length != zLen) { - throw new DimensionMismatchException(d3FdXdYdZ[i][j].length, zLen); - } - - final int jp1 = j + 1; - final double yR = yval[jp1] - yval[j]; - final double xRyR = xR * yR; - for (int k = 0; k < lastK; k++) { - final int kp1 = k + 1; - final double zR = zval[kp1] - zval[k]; - final double xRzR = xR * zR; - final double yRzR = yR * zR; - final double xRyRzR = xR * yRzR; - - final double[] beta = new double[] { - f[i][j][k], f[ip1][j][k], - f[i][jp1][k], f[ip1][jp1][k], - f[i][j][kp1], f[ip1][j][kp1], - f[i][jp1][kp1], f[ip1][jp1][kp1], - - dFdX[i][j][k] * xR, dFdX[ip1][j][k] * xR, - dFdX[i][jp1][k] * xR, dFdX[ip1][jp1][k] * xR, - dFdX[i][j][kp1] * xR, dFdX[ip1][j][kp1] * xR, - dFdX[i][jp1][kp1] * xR, dFdX[ip1][jp1][kp1] * xR, - - dFdY[i][j][k] * yR, dFdY[ip1][j][k] * yR, - dFdY[i][jp1][k] * yR, dFdY[ip1][jp1][k] * yR, - dFdY[i][j][kp1] * yR, dFdY[ip1][j][kp1] * yR, - dFdY[i][jp1][kp1] * yR, dFdY[ip1][jp1][kp1] * yR, - - dFdZ[i][j][k] * zR, dFdZ[ip1][j][k] * zR, - dFdZ[i][jp1][k] * zR, dFdZ[ip1][jp1][k] * zR, - dFdZ[i][j][kp1] * zR, dFdZ[ip1][j][kp1] * zR, - dFdZ[i][jp1][kp1] * zR, dFdZ[ip1][jp1][kp1] * zR, - - d2FdXdY[i][j][k] * xRyR, d2FdXdY[ip1][j][k] * xRyR, - d2FdXdY[i][jp1][k] * xRyR, d2FdXdY[ip1][jp1][k] * xRyR, - d2FdXdY[i][j][kp1] * xRyR, d2FdXdY[ip1][j][kp1] * xRyR, - d2FdXdY[i][jp1][kp1] * xRyR, d2FdXdY[ip1][jp1][kp1] * xRyR, - - d2FdXdZ[i][j][k] * xRzR, d2FdXdZ[ip1][j][k] * xRzR, - d2FdXdZ[i][jp1][k] * xRzR, d2FdXdZ[ip1][jp1][k] * xRzR, - d2FdXdZ[i][j][kp1] * xRzR, d2FdXdZ[ip1][j][kp1] * xRzR, - d2FdXdZ[i][jp1][kp1] * xRzR, d2FdXdZ[ip1][jp1][kp1] * xRzR, - - d2FdYdZ[i][j][k] * yRzR, d2FdYdZ[ip1][j][k] * yRzR, - d2FdYdZ[i][jp1][k] * yRzR, d2FdYdZ[ip1][jp1][k] * yRzR, - d2FdYdZ[i][j][kp1] * yRzR, d2FdYdZ[ip1][j][kp1] * yRzR, - d2FdYdZ[i][jp1][kp1] * yRzR, d2FdYdZ[ip1][jp1][kp1] * yRzR, - - d3FdXdYdZ[i][j][k] * xRyRzR, d3FdXdYdZ[ip1][j][k] * xRyRzR, - d3FdXdYdZ[i][jp1][k] * xRyRzR, d3FdXdYdZ[ip1][jp1][k] * xRyRzR, - d3FdXdYdZ[i][j][kp1] * xRyRzR, d3FdXdYdZ[ip1][j][kp1] * xRyRzR, - d3FdXdYdZ[i][jp1][kp1] * xRyRzR, d3FdXdYdZ[ip1][jp1][kp1] * xRyRzR, - }; - - splines[i][j][k] = new TricubicFunction(computeCoefficients(beta)); - } - } - } - } - - /** - * {@inheritDoc} - * - * @throws OutOfRangeException if any of the variables is outside its interpolation range. - */ - public double value(double x, double y, double z) - throws OutOfRangeException { - final int i = searchIndex(x, xval); - if (i == -1) { - throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]); - } - final int j = searchIndex(y, yval); - if (j == -1) { - throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]); - } - final int k = searchIndex(z, zval); - if (k == -1) { - throw new OutOfRangeException(z, zval[0], zval[zval.length - 1]); - } - - final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); - final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); - final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]); - - return splines[i][j][k].value(xN, yN, zN); - } - - /** - * Indicates whether a point is within the interpolation range. - * - * @param x First coordinate. - * @param y Second coordinate. - * @param z Third coordinate. - * @return {@code true} if (x, y, z) is a valid point. - */ - public boolean isValidPoint(double x, double y, double z) { - if (x < xval[0] || - x > xval[xval.length - 1] || - y < yval[0] || - y > yval[yval.length - 1] || - z < zval[0] || - z > zval[zval.length - 1]) { - return false; - } else { - return true; - } - } - - /** - * @param c Coordinate. - * @param val Coordinate samples. - * @return the index in {@code val} corresponding to the interval containing {@code c}, or {@code -1} - * if {@code c} is out of the range defined by the end values of {@code val}. - */ - private int searchIndex(double c, double[] val) { - if (c < val[0]) { - return -1; - } - - final int max = val.length; - for (int i = 1; i < max; i++) { - if (c <= val[i]) { - return i - 1; - } - } - - return -1; - } - - /** - * Compute the spline coefficients from the list of function values and - * function partial derivatives values at the four corners of a grid - * element. They must be specified in the following order: - * <ul> - * <li>f(0,0,0)</li> - * <li>f(1,0,0)</li> - * <li>f(0,1,0)</li> - * <li>f(1,1,0)</li> - * <li>f(0,0,1)</li> - * <li>f(1,0,1)</li> - * <li>f(0,1,1)</li> - * <li>f(1,1,1)</li> - * - * <li>f<sub>x</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>x</sub>(1,1,1)</li> - * - * <li>f<sub>y</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>y</sub>(1,1,1)</li> - * - * <li>f<sub>z</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>z</sub>(1,1,1)</li> - * - * <li>f<sub>xy</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>xy</sub>(1,1,1)</li> - * - * <li>f<sub>xz</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>xz</sub>(1,1,1)</li> - * - * <li>f<sub>yz</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>yz</sub>(1,1,1)</li> - * - * <li>f<sub>xyz</sub>(0,0,0)</li> - * <li>... <em>(same order as above)</em></li> - * <li>f<sub>xyz</sub>(1,1,1)</li> - * </ul> - * where the subscripts indicate the partial derivative with respect to - * the corresponding variable(s). - * - * @param beta List of function values and function partial derivatives values. - * @return the spline coefficients. - */ - private double[] computeCoefficients(double[] beta) { - final int sz = 64; - final double[] a = new double[sz]; - - for (int i = 0; i < sz; i++) { - double result = 0; - final double[] row = AINV[i]; - for (int j = 0; j < sz; j++) { - result += row[j] * beta[j]; - } - a[i] = result; - } - - return a; - } -} - -/** - * 3D-spline function. - * - */ -class TricubicFunction - implements TrivariateFunction { - /** Number of points. */ - private static final short N = 4; - /** Coefficients */ - private final double[][][] a = new double[N][N][N]; - - /** - * @param aV List of spline coefficients. - */ - public TricubicFunction(double[] aV) { - for (int i = 0; i < N; i++) { - for (int j = 0; j < N; j++) { - for (int k = 0; k < N; k++) { - a[i][j][k] = aV[i + N * (j + N * k)]; - } - } - } - } - - /** - * @param x x-coordinate of the interpolation point. - * @param y y-coordinate of the interpolation point. - * @param z z-coordinate of the interpolation point. - * @return the interpolated value. - * @throws OutOfRangeException if {@code x}, {@code y} or - * {@code z} are not in the interval {@code [0, 1]}. - */ - public double value(double x, double y, double z) - throws OutOfRangeException { - if (x < 0 || x > 1) { - throw new OutOfRangeException(x, 0, 1); - } - if (y < 0 || y > 1) { - throw new OutOfRangeException(y, 0, 1); - } - if (z < 0 || z > 1) { - throw new OutOfRangeException(z, 0, 1); - } - - final double x2 = x * x; - final double x3 = x2 * x; - final double[] pX = { 1, x, x2, x3 }; - - final double y2 = y * y; - final double y3 = y2 * y; - final double[] pY = { 1, y, y2, y3 }; - - final double z2 = z * z; - final double z3 = z2 * z; - final double[] pZ = { 1, z, z2, z3 }; - - double result = 0; - for (int i = 0; i < N; i++) { - for (int j = 0; j < N; j++) { - for (int k = 0; k < N; k++) { - result += a[i][j][k] * pX[i] * pY[j] * pZ[k]; - } - } - } - - return result; - } -} http://git-wip-us.apache.org/repos/asf/commons-math/blob/a7b4803f/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolator.java ---------------------------------------------------------------------- diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolator.java deleted file mode 100644 index 531e736..0000000 --- a/src/main/java/org/apache/commons/math3/analysis/interpolation/TricubicInterpolator.java +++ /dev/null @@ -1,142 +0,0 @@ -/* - * Licensed to the Apache Software Foundation (ASF) under one or more - * contributor license agreements. See the NOTICE file distributed with - * this work for additional information regarding copyright ownership. - * The ASF licenses this file to You under the Apache License, Version 2.0 - * (the "License"); you may not use this file except in compliance with - * the License. You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package org.apache.commons.math3.analysis.interpolation; - -import org.apache.commons.math3.exception.DimensionMismatchException; -import org.apache.commons.math3.exception.NoDataException; -import org.apache.commons.math3.exception.NonMonotonicSequenceException; -import org.apache.commons.math3.exception.NumberIsTooSmallException; -import org.apache.commons.math3.util.MathArrays; - -/** - * Generates a tricubic interpolating function. - * - * @since 3.4 - */ -public class TricubicInterpolator - implements TrivariateGridInterpolator { - /** - * {@inheritDoc} - */ - public TricubicInterpolatingFunction interpolate(final double[] xval, - final double[] yval, - final double[] zval, - final double[][][] fval) - throws NoDataException, NumberIsTooSmallException, - DimensionMismatchException, NonMonotonicSequenceException { - if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) { - throw new NoDataException(); - } - if (xval.length != fval.length) { - throw new DimensionMismatchException(xval.length, fval.length); - } - - MathArrays.checkOrder(xval); - MathArrays.checkOrder(yval); - MathArrays.checkOrder(zval); - - final int xLen = xval.length; - final int yLen = yval.length; - final int zLen = zval.length; - - // Approximation to the partial derivatives using finite differences. - final double[][][] dFdX = new double[xLen][yLen][zLen]; - final double[][][] dFdY = new double[xLen][yLen][zLen]; - final double[][][] dFdZ = new double[xLen][yLen][zLen]; - final double[][][] d2FdXdY = new double[xLen][yLen][zLen]; - final double[][][] d2FdXdZ = new double[xLen][yLen][zLen]; - final double[][][] d2FdYdZ = new double[xLen][yLen][zLen]; - final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen]; - - for (int i = 1; i < xLen - 1; i++) { - if (yval.length != fval[i].length) { - throw new DimensionMismatchException(yval.length, fval[i].length); - } - - final int nI = i + 1; - final int pI = i - 1; - - final double nX = xval[nI]; - final double pX = xval[pI]; - - final double deltaX = nX - pX; - - for (int j = 1; j < yLen - 1; j++) { - if (zval.length != fval[i][j].length) { - throw new DimensionMismatchException(zval.length, fval[i][j].length); - } - - final int nJ = j + 1; - final int pJ = j - 1; - - final double nY = yval[nJ]; - final double pY = yval[pJ]; - - final double deltaY = nY - pY; - final double deltaXY = deltaX * deltaY; - - for (int k = 1; k < zLen - 1; k++) { - final int nK = k + 1; - final int pK = k - 1; - - final double nZ = zval[nK]; - final double pZ = zval[pK]; - - final double deltaZ = nZ - pZ; - - dFdX[i][j][k] = (fval[nI][j][k] - fval[pI][j][k]) / deltaX; - dFdY[i][j][k] = (fval[i][nJ][k] - fval[i][pJ][k]) / deltaY; - dFdZ[i][j][k] = (fval[i][j][nK] - fval[i][j][pK]) / deltaZ; - - final double deltaXZ = deltaX * deltaZ; - final double deltaYZ = deltaY * deltaZ; - - d2FdXdY[i][j][k] = (fval[nI][nJ][k] - fval[nI][pJ][k] - fval[pI][nJ][k] + fval[pI][pJ][k]) / deltaXY; - d2FdXdZ[i][j][k] = (fval[nI][j][nK] - fval[nI][j][pK] - fval[pI][j][nK] + fval[pI][j][pK]) / deltaXZ; - d2FdYdZ[i][j][k] = (fval[i][nJ][nK] - fval[i][nJ][pK] - fval[i][pJ][nK] + fval[i][pJ][pK]) / deltaYZ; - - final double deltaXYZ = deltaXY * deltaZ; - - d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] - - fval[pI][nJ][nK] + fval[pI][pJ][nK] - - fval[nI][nJ][pK] + fval[nI][pJ][pK] + - fval[pI][nJ][pK] - fval[pI][pJ][pK]) / deltaXYZ; - } - } - } - - // Create the interpolating function. - return new TricubicInterpolatingFunction(xval, yval, zval, fval, - dFdX, dFdY, dFdZ, - d2FdXdY, d2FdXdZ, d2FdYdZ, - d3FdXdYdZ) { - @Override - public boolean isValidPoint(double x, double y, double z) { - if (x < xval[1] || - x > xval[xval.length - 2] || - y < yval[1] || - y > yval[yval.length - 2] || - z < zval[1] || - z > zval[zval.length - 2]) { - return false; - } else { - return true; - } - } - }; - } -}
