-
Notifications
You must be signed in to change notification settings - Fork 41
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Added bilinear interpolation for 2D grids.
Github: fixes #48.
- Loading branch information
Showing
7 changed files
with
763 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
138 changes: 138 additions & 0 deletions
138
...re/src/main/java/org/hipparchus/analysis/interpolation/BilinearInterpolatingFunction.java
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,138 @@ | ||
| /* | ||
| * Licensed to the Hipparchus project under one or more | ||
| * contributor license agreements. See the NOTICE file distributed with | ||
| * this work for additional information regarding copyright ownership. | ||
| * The Hipparchus project 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.hipparchus.analysis.interpolation; | ||
|
|
||
| import java.io.Serializable; | ||
|
|
||
| import org.hipparchus.analysis.BivariateFunction; | ||
| import org.hipparchus.exception.MathIllegalArgumentException; | ||
|
|
||
| /** | ||
| * Interpolate grid data using bi-linear interpolation. | ||
| * <p> | ||
| * This interpolator is thread-safe. | ||
| * </p> | ||
| * @since 1.4 | ||
| */ | ||
| public class BilinearInterpolatingFunction implements BivariateFunction, Serializable { | ||
|
|
||
| /** Serializable UID. */ | ||
| private static final long serialVersionUID = 20180926L; | ||
|
|
||
| /** Grid along the x axis. */ | ||
| private final GridAxis xGrid; | ||
|
|
||
| /** Grid size along the x axis. */ | ||
| private final int xSize; | ||
|
|
||
| /** Grid along the y axis. */ | ||
| private final GridAxis yGrid; | ||
|
|
||
| /** Grid size along the y axis. */ | ||
| private final int ySize; | ||
|
|
||
| /** Values of the interpolation points on all the grid knots (in a flatten array). */ | ||
| private final double[] fVal; | ||
|
|
||
| /** Simple constructor. | ||
| * @param xVal All the x-coordinates of the interpolation points, sorted | ||
| * in increasing order. | ||
| * @param yVal All the y-coordinates of the interpolation points, sorted | ||
| * in increasing order. | ||
| * @param fVal The values of the interpolation points on all the grid knots: | ||
| * {@code fVal[i][j] = f(xVal[i], yVal[j])}. | ||
| * @exception MathIllegalArgumentException if grid size is smaller than 2 | ||
| * or if the grid is not sorted in strict increasing order | ||
| */ | ||
| public BilinearInterpolatingFunction(final double[] xVal, final double[] yVal, | ||
| final double[][] fVal) | ||
| throws MathIllegalArgumentException { | ||
| this.xGrid = new GridAxis(xVal, 2); | ||
| this.xSize = xVal.length; | ||
| this.yGrid = new GridAxis(yVal, 2); | ||
| this.ySize = yVal.length; | ||
| this.fVal = new double[xSize * ySize]; | ||
| int k = 0; | ||
| for (int i = 0; i < xSize; ++i) { | ||
| final double[] fi = fVal[i]; | ||
| for (int j = 0; j < ySize; ++j) { | ||
| this.fVal[k++] = fi[j]; | ||
| } | ||
| } | ||
| } | ||
|
|
||
| /** Get the lowest grid x coordinate. | ||
| * @return lowest grid x coordinate | ||
| */ | ||
| public double getXInf() { | ||
| return xGrid.node(0); | ||
| } | ||
|
|
||
| /** Get the highest grid x coordinate. | ||
| * @return highest grid x coordinate | ||
| */ | ||
| public double getXSup() { | ||
| return xGrid.node(xGrid.size() - 1); | ||
| } | ||
|
|
||
| /** Get the lowest grid y coordinate. | ||
| * @return lowest grid y coordinate | ||
| */ | ||
| public double getYInf() { | ||
| return yGrid.node(0); | ||
| } | ||
|
|
||
| /** Get the highest grid y coordinate. | ||
| * @return highest grid y coordinate | ||
| */ | ||
| public double getYSup() { | ||
| return yGrid.node(yGrid.size() - 1); | ||
| } | ||
|
|
||
| /** {@inheritDoc} */ | ||
| @Override | ||
| public double value(final double x, final double y) { | ||
|
|
||
| // get the interpolation nodes | ||
| final int i = xGrid.interpolationIndex(x); | ||
| final int j = yGrid.interpolationIndex(y); | ||
| final double x0 = xGrid.node(i); | ||
| final double x1 = xGrid.node(i + 1); | ||
| final double y0 = yGrid.node(j); | ||
| final double y1 = yGrid.node(j + 1); | ||
|
|
||
| // get the function values at interpolation nodes | ||
| final int k0 = i * ySize + j; | ||
| final int k1 = k0 + ySize; | ||
| final double z00 = fVal[k0]; | ||
| final double z01 = fVal[k0 + 1]; | ||
| final double z10 = fVal[k1]; | ||
| final double z11 = fVal[k1 + 1]; | ||
|
|
||
| // interpolate | ||
| final double dx0 = x - x0; | ||
| final double dx1 = x1 - x; | ||
| final double dx10 = x1 - x0; | ||
| final double dy0 = y - y0; | ||
| final double dy1 = y1 - y; | ||
| final double dy10 = y1 - y0; | ||
| return (dx0 * (dy0 * z11 + dy1 * z10) + dx1 * (dy0 * z01 + dy1 * z00)) / | ||
| (dx10 * dy10); | ||
|
|
||
| } | ||
|
|
||
| } |
35 changes: 35 additions & 0 deletions
35
...archus-core/src/main/java/org/hipparchus/analysis/interpolation/BilinearInterpolator.java
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,35 @@ | ||
| /* | ||
| * Licensed to the Hipparchus project under one or more | ||
| * contributor license agreements. See the NOTICE file distributed with | ||
| * this work for additional information regarding copyright ownership. | ||
| * The Hipparchus project 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.hipparchus.analysis.interpolation; | ||
|
|
||
| import org.hipparchus.exception.MathIllegalArgumentException; | ||
|
|
||
| /** | ||
| * Interpolate grid data using bi-linear interpolation. | ||
| * @since 1.4 | ||
| */ | ||
| public class BilinearInterpolator implements BivariateGridInterpolator { | ||
|
|
||
| /** {@inheritDoc} */ | ||
| @Override | ||
| public BilinearInterpolatingFunction interpolate(final double[] xval, final double[] yval, | ||
| final double[][] fval) | ||
| throws MathIllegalArgumentException { | ||
| return new BilinearInterpolatingFunction(xval, yval, fval); | ||
| } | ||
|
|
||
| } |
173 changes: 173 additions & 0 deletions
173
hipparchus-core/src/main/java/org/hipparchus/analysis/interpolation/GridAxis.java
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,173 @@ | ||
| /* | ||
| * Licensed to the Hipparchus project under one or more | ||
| * contributor license agreements. See the NOTICE file distributed with | ||
| * this work for additional information regarding copyright ownership. | ||
| * The Hipparchus project 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.hipparchus.analysis.interpolation; | ||
|
|
||
| import java.io.Serializable; | ||
| import java.util.concurrent.atomic.AtomicInteger; | ||
|
|
||
| import org.hipparchus.exception.LocalizedCoreFormats; | ||
| import org.hipparchus.exception.MathIllegalArgumentException; | ||
| import org.hipparchus.util.FastMath; | ||
| import org.hipparchus.util.MathArrays; | ||
|
|
||
| /** | ||
| * Helper for finding interpolation nodes along one axis of grid data. | ||
| * <p> | ||
| * This class is intended to be used for interpolating inside grids. | ||
| * It works on any sorted data without duplication and size at least | ||
| * {@code n} where {@code n} is the number of points required for | ||
| * interpolation (i.e. 2 for linear interpolation, 3 for quadratic...) | ||
| * <p> | ||
| * </p> | ||
| * The method uses linear interpolation to select the nodes indices. | ||
| * It should be O(1) for sufficiently regular data, therefore much faster | ||
| * than bisection. It also features caching, which improves speed when | ||
| * interpolating several points in raw in the close locations, i.e. when | ||
| * successive calls have a high probability to return the same interpolation | ||
| * nodes. This occurs for example when scanning with small steps a loose | ||
| * grid. The method also works on non-regular grids, but may be slower in | ||
| * this case. | ||
| * </p> | ||
| * <p> | ||
| * This class is thread-safe. | ||
| * </p> | ||
| * @since 1.4 | ||
| */ | ||
| public class GridAxis implements Serializable { | ||
|
|
||
| /** Serializable UID. */ | ||
| private static final long serialVersionUID = 20180926L; | ||
|
|
||
| /** All the coordinates of the interpolation points, sorted in increasing order. */ | ||
| private final double[] grid; | ||
|
|
||
| /** Number of points required for interpolation. */ | ||
| private int n; | ||
|
|
||
| /** Cached value of last x index. */ | ||
| private final AtomicInteger cache; | ||
|
|
||
| /** Simple constructor. | ||
| * @param grid coordinates of the interpolation points, sorted in increasing order | ||
| * @param n number of points required for interpolation, i.e. 2 for linear, 3 | ||
| * for quadratic... | ||
| * @exception MathIllegalArgumentException if grid size is smaller than {@code n} | ||
| * or if the grid is not sorted in strict increasing order | ||
| */ | ||
| public GridAxis(final double[] grid, final int n) | ||
| throws MathIllegalArgumentException { | ||
|
|
||
| // safety checks | ||
| if (grid.length < n) { | ||
| throw new MathIllegalArgumentException(LocalizedCoreFormats.INSUFFICIENT_DIMENSION, | ||
| grid.length, n); | ||
| } | ||
| MathArrays.checkOrder(grid); | ||
|
|
||
| this.grid = grid.clone(); | ||
| this.n = n; | ||
| this.cache = new AtomicInteger(0); | ||
|
|
||
| } | ||
|
|
||
| /** Get the number of points of the grid. | ||
| * @return number of points of the grid | ||
| */ | ||
| public int size() { | ||
| return grid.length; | ||
| } | ||
|
|
||
| /** Get the number of points required for interpolation. | ||
| * @return number of points required for interpolation | ||
| */ | ||
| public int getN() { | ||
| return n; | ||
| } | ||
|
|
||
| /** Get the interpolation node at specified index. | ||
| * @param index node index | ||
| * @return coordinate of the node at specified index | ||
| */ | ||
| public double node(final int index) { | ||
| return grid[index]; | ||
| } | ||
|
|
||
| /** Get the index of the first interpolation node for some coordinate along the grid. | ||
| * <p> | ||
| * The index return is the one for the lowest interpolation node suitabme for | ||
| * {@code t}. This means that if {@code i} is returned the nodes to use for | ||
| * interpolation at coordinate {@code t} are at indices {@code i}, {@code i+1}, | ||
| * ..., {@code i+n-1}, where {@code n} is the number of points required for | ||
| * interpolation passed at construction. | ||
| * </p> | ||
| * <p> | ||
| * The index is selected in order to have the array as balanced as possible | ||
| * around {@code t}. This means that if {@code n} if even and {@code t} is | ||
| * sufficiently far from the grid endpoints, the nodes at indices {@code i} to | ||
| * {code i+n/2} will be smaller than {@code t} and nodes at indices {@code i+n/2+1} | ||
| * to {@code i+n-1} will be larger than {@code t}. If {@code t} is close to endpoints | ||
| * but still within the grid, the returned node will be unbalanced. If {@code t} | ||
| * is outside of the grid, all nodes will be on the same side with respect to | ||
| * {@code t}. It is <em>not</em> an error to call this method with {@code t} outside | ||
| * of the grid, it simply implies that the interpolation will become an extrapolation | ||
| * and accuracy will decrease as {@code t} goes farther from the grid points. This | ||
| * is intended so interpolation does not fail near the end of the grid. | ||
| * </p> | ||
| * @param t coordinate of the point to interpolate | ||
| * @return index {@code i} such {@link #node(int) node(i)}, {@link #node(int) node(i+1)}, | ||
| * ... {@link #node(int) node(i+n-1)} can be used for interpolating a value at | ||
| * coordinate {@code t} | ||
| * @since 1.4 | ||
| */ | ||
| public int interpolationIndex(final double t) { | ||
|
|
||
| final int middleOffset = (n - 1) / 2; | ||
|
|
||
| // first try to simply reuse the last used index, | ||
| // for faster return in a common case | ||
| final int cached = cache.get(); | ||
| final int middle = cached + middleOffset; | ||
| if (grid[middle] <= t && t < grid[middle + 1]) { | ||
| return cached; | ||
| } | ||
|
|
||
| // we need to find a new index | ||
| int iInf = middleOffset; | ||
| double aInf = grid[iInf]; | ||
| int iSup = grid.length - (n - 1) + middleOffset; | ||
| double aSup = grid[iSup]; | ||
| while (iSup - iInf > 1) { | ||
| final int iInterp = (int) ((iInf * (aSup - t) + iSup * (t - aInf)) / (aSup - aInf)); | ||
| final int iMed = FastMath.max(iInf + 1, FastMath.min(iInterp, iSup - 1)); | ||
| if (t < grid[iMed]) { | ||
| // keeps looking in the lower part of the grid | ||
| iSup = iMed; | ||
| aSup = grid[iSup]; | ||
| } else { | ||
| // keeps looking in the upper part of the grid | ||
| iInf = iMed; | ||
| aInf = grid[iInf]; | ||
| } | ||
| } | ||
|
|
||
| final int newCached = iInf - middleOffset; | ||
| cache.compareAndSet(cached, newCached); | ||
| return newCached; | ||
|
|
||
| } | ||
|
|
||
| } |
Oops, something went wrong.