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.analysis.TrivariateRealFunction;
020 import org.apache.commons.math.exception.DimensionMismatchException;
021 import org.apache.commons.math.exception.NoDataException;
022 import org.apache.commons.math.exception.OutOfRangeException;
023 import org.apache.commons.math.util.MathUtils;
024
025 /**
026 * Function that implements the
027 * <a href="http://en.wikipedia.org/wiki/Tricubic_interpolation">
028 * tricubic spline interpolation</a>, as proposed in
029 * <quote>
030 * Tricubic interpolation in three dimensions<br/>
031 * F. Lekien and J. Marsden<br/>
032 * <em>Int. J. Numer. Meth. Engng</em> 2005; <b>63</b>:455-471
033 * </quote>
034 *
035 * @version $Revision$ $Date$
036 * @since 2.2
037 */
038 public class TricubicSplineInterpolatingFunction
039 implements TrivariateRealFunction {
040 /**
041 * Matrix to compute the spline coefficients from the function values
042 * and function derivatives values
043 */
044 private static final double[][] AINV = {
045 { 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 },
046 { 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 },
047 { -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 },
048 { 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 },
049 { 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 },
050 { 0,0,0,0,0,0,0,0,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 },
051 { 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 },
052 { 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 },
053 { -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 },
054 { 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 },
055 { 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 },
056 { -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 },
057 { 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 },
058 { 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 },
059 { -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 },
060 { 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 },
061 { 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 },
062 { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
063 { 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 },
064 { 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 },
065 { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
066 { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
067 { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
068 { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
069 { 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 },
070 { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
071 { 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 },
072 { 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 },
073 { 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 },
074 { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,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 },
075 { 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 },
076 { 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 },
077 {-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 },
078 { 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 },
079 { 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 },
080 { -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 },
081 { 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 },
082 { 0,0,0,0,0,0,0,0,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 },
083 { 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 },
084 { 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 },
085 { 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 },
086 { 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 },
087 { -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 },
088 { 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 },
089 { -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 },
090 { 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 },
091 { 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 },
092 { -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 },
093 { 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 },
094 { 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 },
095 { -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 },
096 { 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 },
097 { 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 },
098 { 0,0,0,0,0,0,0,0,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 },
099 { 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 },
100 { 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 },
101 { -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 },
102 { 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 },
103 { 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 },
104 { -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 },
105 { 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 },
106 { 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 },
107 { -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 },
108 { 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 }
109 };
110
111 /** Samples x-coordinates */
112 private final double[] xval;
113 /** Samples y-coordinates */
114 private final double[] yval;
115 /** Samples z-coordinates */
116 private final double[] zval;
117 /** Set of cubic splines pacthing the whole data grid */
118 private final TricubicSplineFunction[][][] splines;
119
120 /**
121 * @param x Sample values of the x-coordinate, in increasing order.
122 * @param y Sample values of the y-coordinate, in increasing order.
123 * @param z Sample values of the y-coordinate, in increasing order.
124 * @param f Values of the function on every grid point.
125 * @param dFdX Values of the partial derivative of function with respect
126 * to x on every grid point.
127 * @param dFdY Values of the partial derivative of function with respect
128 * to y on every grid point.
129 * @param dFdZ Values of the partial derivative of function with respect
130 * to z on every grid point.
131 * @param d2FdXdY Values of the cross partial derivative of function on
132 * every grid point.
133 * @param d2FdXdZ Values of the cross partial derivative of function on
134 * every grid point.
135 * @param d2FdYdZ Values of the cross partial derivative of function on
136 * every grid point.
137 * @param d3FdXdYdZ Values of the cross partial derivative of function on
138 * every grid point.
139 * @throws NoDataException if any of the arrays has zero length.
140 * @throws DimensionMismatchException if the various arrays do not contain
141 * the expected number of elements.
142 * @throws IllegalArgumentException if {@code x}, {@code y} or {@code z}
143 * are not strictly increasing.
144 */
145 public TricubicSplineInterpolatingFunction(double[] x,
146 double[] y,
147 double[] z,
148 double[][][] f,
149 double[][][] dFdX,
150 double[][][] dFdY,
151 double[][][] dFdZ,
152 double[][][] d2FdXdY,
153 double[][][] d2FdXdZ,
154 double[][][] d2FdYdZ,
155 double[][][] d3FdXdYdZ) {
156 final int xLen = x.length;
157 final int yLen = y.length;
158 final int zLen = z.length;
159
160 if (xLen == 0 || yLen == 0 || z.length == 0 || f.length == 0 || f[0].length == 0) {
161 throw new NoDataException();
162 }
163 if (xLen != f.length) {
164 throw new DimensionMismatchException(xLen, f.length);
165 }
166 if (xLen != dFdX.length) {
167 throw new DimensionMismatchException(xLen, dFdX.length);
168 }
169 if (xLen != dFdY.length) {
170 throw new DimensionMismatchException(xLen, dFdY.length);
171 }
172 if (xLen != dFdZ.length) {
173 throw new DimensionMismatchException(xLen, dFdZ.length);
174 }
175 if (xLen != d2FdXdY.length) {
176 throw new DimensionMismatchException(xLen, d2FdXdY.length);
177 }
178 if (xLen != d2FdXdZ.length) {
179 throw new DimensionMismatchException(xLen, d2FdXdZ.length);
180 }
181 if (xLen != d2FdYdZ.length) {
182 throw new DimensionMismatchException(xLen, d2FdYdZ.length);
183 }
184 if (xLen != d3FdXdYdZ.length) {
185 throw new DimensionMismatchException(xLen, d3FdXdYdZ.length);
186 }
187
188 MathUtils.checkOrder(x);
189 MathUtils.checkOrder(y);
190 MathUtils.checkOrder(z);
191
192 xval = x.clone();
193 yval = y.clone();
194 zval = z.clone();
195
196 final int lastI = xLen - 1;
197 final int lastJ = yLen - 1;
198 final int lastK = zLen - 1;
199 splines = new TricubicSplineFunction[lastI][lastJ][lastK];
200
201 for (int i = 0; i < lastI; i++) {
202 if (f[i].length != yLen) {
203 throw new DimensionMismatchException(f[i].length, yLen);
204 }
205 if (dFdX[i].length != yLen) {
206 throw new DimensionMismatchException(dFdX[i].length, yLen);
207 }
208 if (dFdY[i].length != yLen) {
209 throw new DimensionMismatchException(dFdY[i].length, yLen);
210 }
211 if (dFdZ[i].length != yLen) {
212 throw new DimensionMismatchException(dFdZ[i].length, yLen);
213 }
214 if (d2FdXdY[i].length != yLen) {
215 throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
216 }
217 if (d2FdXdZ[i].length != yLen) {
218 throw new DimensionMismatchException(d2FdXdZ[i].length, yLen);
219 }
220 if (d2FdYdZ[i].length != yLen) {
221 throw new DimensionMismatchException(d2FdYdZ[i].length, yLen);
222 }
223 if (d3FdXdYdZ[i].length != yLen) {
224 throw new DimensionMismatchException(d3FdXdYdZ[i].length, yLen);
225 }
226
227 final int ip1 = i + 1;
228 for (int j = 0; j < lastJ; j++) {
229 if (f[i][j].length != zLen) {
230 throw new DimensionMismatchException(f[i][j].length, zLen);
231 }
232 if (dFdX[i][j].length != zLen) {
233 throw new DimensionMismatchException(dFdX[i][j].length, zLen);
234 }
235 if (dFdY[i][j].length != zLen) {
236 throw new DimensionMismatchException(dFdY[i][j].length, zLen);
237 }
238 if (dFdZ[i][j].length != zLen) {
239 throw new DimensionMismatchException(dFdZ[i][j].length, zLen);
240 }
241 if (d2FdXdY[i][j].length != zLen) {
242 throw new DimensionMismatchException(d2FdXdY[i][j].length, zLen);
243 }
244 if (d2FdXdZ[i][j].length != zLen) {
245 throw new DimensionMismatchException(d2FdXdZ[i][j].length, zLen);
246 }
247 if (d2FdYdZ[i][j].length != zLen) {
248 throw new DimensionMismatchException(d2FdYdZ[i][j].length, zLen);
249 }
250 if (d3FdXdYdZ[i][j].length != zLen) {
251 throw new DimensionMismatchException(d3FdXdYdZ[i][j].length, zLen);
252 }
253
254 final int jp1 = j + 1;
255 for (int k = 0; k < lastK; k++) {
256 final int kp1 = k + 1;
257
258 final double[] beta = new double[] {
259 f[i][j][k], f[ip1][j][k],
260 f[i][jp1][k], f[ip1][jp1][k],
261 f[i][j][kp1], f[ip1][j][kp1],
262 f[i][jp1][kp1], f[ip1][jp1][kp1],
263
264 dFdX[i][j][k], dFdX[ip1][j][k],
265 dFdX[i][jp1][k], dFdX[ip1][jp1][k],
266 dFdX[i][j][kp1], dFdX[ip1][j][kp1],
267 dFdX[i][jp1][kp1], dFdX[ip1][jp1][kp1],
268
269 dFdY[i][j][k], dFdY[ip1][j][k],
270 dFdY[i][jp1][k], dFdY[ip1][jp1][k],
271 dFdY[i][j][kp1], dFdY[ip1][j][kp1],
272 dFdY[i][jp1][kp1], dFdY[ip1][jp1][kp1],
273
274 dFdZ[i][j][k], dFdZ[ip1][j][k],
275 dFdZ[i][jp1][k], dFdZ[ip1][jp1][k],
276 dFdZ[i][j][kp1], dFdZ[ip1][j][kp1],
277 dFdZ[i][jp1][kp1], dFdZ[ip1][jp1][kp1],
278
279 d2FdXdY[i][j][k], d2FdXdY[ip1][j][k],
280 d2FdXdY[i][jp1][k], d2FdXdY[ip1][jp1][k],
281 d2FdXdY[i][j][kp1], d2FdXdY[ip1][j][kp1],
282 d2FdXdY[i][jp1][kp1], d2FdXdY[ip1][jp1][kp1],
283
284 d2FdXdZ[i][j][k], d2FdXdZ[ip1][j][k],
285 d2FdXdZ[i][jp1][k], d2FdXdZ[ip1][jp1][k],
286 d2FdXdZ[i][j][kp1], d2FdXdZ[ip1][j][kp1],
287 d2FdXdZ[i][jp1][kp1], d2FdXdZ[ip1][jp1][kp1],
288
289 d2FdYdZ[i][j][k], d2FdYdZ[ip1][j][k],
290 d2FdYdZ[i][jp1][k], d2FdYdZ[ip1][jp1][k],
291 d2FdYdZ[i][j][kp1], d2FdYdZ[ip1][j][kp1],
292 d2FdYdZ[i][jp1][kp1], d2FdYdZ[ip1][jp1][kp1],
293
294 d3FdXdYdZ[i][j][k], d3FdXdYdZ[ip1][j][k],
295 d3FdXdYdZ[i][jp1][k], d3FdXdYdZ[ip1][jp1][k],
296 d3FdXdYdZ[i][j][kp1], d3FdXdYdZ[ip1][j][kp1],
297 d3FdXdYdZ[i][jp1][kp1], d3FdXdYdZ[ip1][jp1][kp1],
298 };
299
300 splines[i][j][k] = new TricubicSplineFunction(computeSplineCoefficients(beta));
301 }
302 }
303 }
304 }
305
306 /**
307 * {@inheritDoc}
308 */
309 public double value(double x, double y, double z) {
310 final int i = searchIndex(x, xval);
311 if (i == -1) {
312 throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]);
313 }
314 final int j = searchIndex(y, yval);
315 if (j == -1) {
316 throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]);
317 }
318 final int k = searchIndex(z, zval);
319 if (k == -1) {
320 throw new OutOfRangeException(z, zval[0], zval[zval.length - 1]);
321 }
322
323 final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
324 final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
325 final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]);
326
327 return splines[i][j][k].value(xN, yN, zN);
328 }
329
330 /**
331 * @param c Coordinate.
332 * @param val Coordinate samples.
333 * @return the index in {@code val} corresponding to the interval
334 * containing {@code c}, or {@code -1} if {@code c} is out of the
335 * range defined by the end values of {@code val}.
336 */
337 private int searchIndex(double c, double[] val) {
338 if (c < val[0]) {
339 return -1;
340 }
341
342 final int max = val.length;
343 for (int i = 1; i < max; i++) {
344 if (c <= val[i]) {
345 return i - 1;
346 }
347 }
348
349 return -1;
350 }
351
352 /**
353 * Compute the spline coefficients from the list of function values and
354 * function partial derivatives values at the four corners of a grid
355 * element. They must be specified in the following order:
356 * <ul>
357 * <li>f(0,0,0)</li>
358 * <li>f(1,0,0)</li>
359 * <li>f(0,1,0)</li>
360 * <li>f(1,1,0)</li>
361 * <li>f(0,0,1)</li>
362 * <li>f(1,0,1)</li>
363 * <li>f(0,1,1)</li>
364 * <li>f(1,1,1)</li>
365 *
366 * <li>f<sub>x</sub>(0,0,0)</li>
367 * <li>... <em>(same order as above)</em></li>
368 * <li>f<sub>x</sub>(1,1,1)</li>
369 *
370 * <li>f<sub>y</sub>(0,0,0)</li>
371 * <li>... <em>(same order as above)</em></li>
372 * <li>f<sub>y</sub>(1,1,1)</li>
373 *
374 * <li>f<sub>z</sub>(0,0,0)</li>
375 * <li>... <em>(same order as above)</em></li>
376 * <li>f<sub>z</sub>(1,1,1)</li>
377 *
378 * <li>f<sub>xy</sub>(0,0,0)</li>
379 * <li>... <em>(same order as above)</em></li>
380 * <li>f<sub>xy</sub>(1,1,1)</li>
381 *
382 * <li>f<sub>xz</sub>(0,0,0)</li>
383 * <li>... <em>(same order as above)</em></li>
384 * <li>f<sub>xz</sub>(1,1,1)</li>
385 *
386 * <li>f<sub>yz</sub>(0,0,0)</li>
387 * <li>... <em>(same order as above)</em></li>
388 * <li>f<sub>yz</sub>(1,1,1)</li>
389 *
390 * <li>f<sub>xyz</sub>(0,0,0)</li>
391 * <li>... <em>(same order as above)</em></li>
392 * <li>f<sub>xyz</sub>(1,1,1)</li>
393 * </ul>
394 * where the subscripts indicate the partial derivative with respect to
395 * the corresponding variable(s).
396 *
397 * @param beta List of function values and function partial derivatives
398 * values.
399 * @return the spline coefficients.
400 */
401 private double[] computeSplineCoefficients(double[] beta) {
402 final int sz = 64;
403 final double[] a = new double[sz];
404
405 for (int i = 0; i < sz; i++) {
406 double result = 0;
407 final double[] row = AINV[i];
408 for (int j = 0; j < sz; j++) {
409 result += row[j] * beta[j];
410 }
411 a[i] = result;
412 }
413
414 return a;
415 }
416 }
417
418 /**
419 * 3D-spline function.
420 *
421 * @version $Revision$ $Date$
422 */
423 class TricubicSplineFunction
424 implements TrivariateRealFunction {
425 /** Number of points. */
426 private static final short N = 4;
427 /** Coefficients */
428 private final double[][][] a = new double[N][N][N];
429
430 /**
431 * @param aV List of spline coefficients.
432 */
433 public TricubicSplineFunction(double[] aV) {
434 for (int i = 0; i < N; i++) {
435 for (int j = 0; j < N; j++) {
436 for (int k = 0; k < N; k++) {
437 a[i][j][k] = aV[i + N * (j + N * k)];
438 }
439 }
440 }
441 }
442
443 /**
444 * @param x x-coordinate of the interpolation point.
445 * @param y y-coordinate of the interpolation point.
446 * @param z z-coordinate of the interpolation point.
447 * @return the interpolated value.
448 */
449 public double value(double x, double y, double z) {
450 if (x < 0 || x > 1) {
451 throw new OutOfRangeException(x, 0, 1);
452 }
453 if (y < 0 || y > 1) {
454 throw new OutOfRangeException(y, 0, 1);
455 }
456 if (z < 0 || z > 1) {
457 throw new OutOfRangeException(z, 0, 1);
458 }
459
460 final double x2 = x * x;
461 final double x3 = x2 * x;
462 final double[] pX = { 1, x, x2, x3 };
463
464 final double y2 = y * y;
465 final double y3 = y2 * y;
466 final double[] pY = { 1, y, y2, y3 };
467
468 final double z2 = z * z;
469 final double z3 = z2 * z;
470 final double[] pZ = { 1, z, z2, z3 };
471
472 double result = 0;
473 for (int i = 0; i < N; i++) {
474 for (int j = 0; j < N; j++) {
475 for (int k = 0; k < N; k++) {
476 result += a[i][j][k] * pX[i] * pY[j] * pZ[k];
477 }
478 }
479 }
480
481 return result;
482 }
483 }