From 5c3a0448de36013cceef0a1fee0cef5bf45383ff Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Sat, 15 Mar 2014 11:05:01 +0100
Subject: [PATCH] Added a tolerance around tile for elevation interpolation.

Elevation is going to be interpolated at tiles entry and exit points,
which are computed from geodetic conversions. The result points may be
very slightly out of tiles, at numerical accuracy level. These points
should nevertheless be allowed to be interpolated using the closest
pixel data, so a tolerance was needed.

The SimpleTile class uses a fixed tolerance, arbitrarily set to 1/8
pixel.
---
 .../orekit/rugged/core/dem/SimpleTile.java    | 25 ++++--
 .../java/org/orekit/rugged/core/dem/Tile.java | 10 ++-
 .../rugged/core/dem/SimpleTileTest.java       | 85 ++++++++++++++++---
 3 files changed, 99 insertions(+), 21 deletions(-)

diff --git a/rugged-core/src/main/java/org/orekit/rugged/core/dem/SimpleTile.java b/rugged-core/src/main/java/org/orekit/rugged/core/dem/SimpleTile.java
index 9503e0b1..832d8ae2 100644
--- a/rugged-core/src/main/java/org/orekit/rugged/core/dem/SimpleTile.java
+++ b/rugged-core/src/main/java/org/orekit/rugged/core/dem/SimpleTile.java
@@ -27,6 +27,9 @@ import org.orekit.rugged.api.RuggedMessages;
  */
 public class SimpleTile implements Tile {
 
+    /** Tolerance used to interpolate points slightly out of tile (in pixels). */
+    private static final double TOLERANCE = 1.0 / 8.0;
+
     /** Minimum latitude. */
     private double minLatitude;
 
@@ -183,17 +186,20 @@ public class SimpleTile implements Tile {
         return elevations[latitudeIndex * getLongitudeColumns() + longitudeIndex];
     }
 
-    /** {@inheritDoc} */
+    /** {@inheritDoc}
+     * <p>
+     * This classes uses an arbitrary 1/8 pixel tolerance for interpolating
+     * slightly out of tile points.
+     * </p>
+     */
     @Override
-    public double interpolateElevation(double latitude, double longitude)
+    public double interpolateElevation(final double latitude, final double longitude)
         throws RuggedException {
 
         final double doubleLatitudeIndex  = getDoubleLatitudeIndex(latitude);
         final double doubleLongitudeIndex = getDoubleLontitudeIndex(longitude);
-        final int    latitudeIndex        = (int) FastMath.floor(doubleLatitudeIndex);
-        final int    longitudeIndex       = (int) FastMath.floor(doubleLongitudeIndex);
-        if (latitudeIndex  < 0 || latitudeIndex  >= (latitudeRows - 2) ||
-            longitudeIndex < 0 || longitudeIndex >= (longitudeColumns - 2)) {
+        if (doubleLatitudeIndex  < -TOLERANCE || doubleLatitudeIndex  >= (latitudeRows - 1 + TOLERANCE) ||
+            doubleLongitudeIndex < -TOLERANCE || doubleLongitudeIndex >= (longitudeColumns - 1 + TOLERANCE)) {
             throw new RuggedException(RuggedMessages.OUT_OF_TILE_ANGLES,
                                       FastMath.toDegrees(latitude),
                                       FastMath.toDegrees(longitude),
@@ -203,6 +209,13 @@ public class SimpleTile implements Tile {
                                       FastMath.toDegrees(getMaximumLongitude()));
         }
 
+        final int latitudeIndex  = FastMath.max(0,
+                                                FastMath.min(latitudeRows - 2,
+                                                             (int) FastMath.floor(doubleLatitudeIndex)));
+        final int longitudeIndex = FastMath.max(0,
+                                                FastMath.min(longitudeColumns - 2,
+                                                             (int) FastMath.floor(doubleLongitudeIndex)));
+
         // bi-linear interpolation
         final double dLat = doubleLatitudeIndex  - latitudeIndex;
         final double dLon = doubleLongitudeIndex - longitudeIndex;
diff --git a/rugged-core/src/main/java/org/orekit/rugged/core/dem/Tile.java b/rugged-core/src/main/java/org/orekit/rugged/core/dem/Tile.java
index 0cc78fac..1e02b64f 100644
--- a/rugged-core/src/main/java/org/orekit/rugged/core/dem/Tile.java
+++ b/rugged-core/src/main/java/org/orekit/rugged/core/dem/Tile.java
@@ -115,10 +115,18 @@ public interface Tile extends UpdatableTile {
         throws RuggedException;
 
     /** Interpolate elevation.
+     * <p>
+     * In order to cope with numerical accuracy issues when computing
+     * points at tile boundary, a slight tolerance (typically 1/8 pixel)
+     * around the tile is allowed. Elevation can therefore be interpolated
+     * (really extrapolated in this case) even for points slightly overshooting
+     * tile boundaries, using the closest tile pixel. Attempting to interpolate
+     * too far from the tile will trigger an exception.
+     * </p>
      * @param latitude ground point latitude
      * @param longitude ground point longitude
      * @return interpolated elevation
-     * @exception RuggedException if point does not lie within the tile
+     * @exception RuggedException if point is farthest from the tile than the tolerance
      */
     double interpolateElevation(double latitude, double longitude)
         throws RuggedException;
diff --git a/rugged-core/src/test/java/org/orekit/rugged/core/dem/SimpleTileTest.java b/rugged-core/src/test/java/org/orekit/rugged/core/dem/SimpleTileTest.java
index 2e5844f1..bab7daf5 100644
--- a/rugged-core/src/test/java/org/orekit/rugged/core/dem/SimpleTileTest.java
+++ b/rugged-core/src/test/java/org/orekit/rugged/core/dem/SimpleTileTest.java
@@ -16,6 +16,7 @@
  */
 package org.orekit.rugged.core.dem;
 
+import org.apache.commons.math3.util.FastMath;
 import org.junit.Assert;
 import org.junit.Test;
 import org.orekit.rugged.api.RuggedException;
@@ -95,20 +96,7 @@ public class SimpleTileTest {
     }
 
     @Test
-    public void testInterpolation() throws RuggedException {
-        SimpleTile tile = new SimpleTileFactory().createTile();
-        tile.setGeometry(0.0, 0.0, 1.0, 1.0, 50, 50);
-        tile.setElevation(20, 14,  91.0);
-        tile.setElevation(20, 15, 210.0);
-        tile.setElevation(21, 14, 162.0);
-        tile.setElevation(21, 15,  95.0);
-        Assert.assertEquals(150.5, tile.interpolateElevation(20.0, 14.5), 1.0e-10);
-        Assert.assertEquals(128.5, tile.interpolateElevation(21.0, 14.5), 1.0e-10);
-        Assert.assertEquals(146.1, tile.interpolateElevation(20.2, 14.5), 1.0e-10);
-    }
-
-    @Test
-    public void testOutOfBounds() throws RuggedException {
+    public void testOutOfBoundsIndices() throws RuggedException {
 
         SimpleTile tile = new SimpleTileFactory().createTile();
         tile.setGeometry(1.0, 2.0, 0.1, 0.2, 100, 200);
@@ -131,4 +119,73 @@ public class SimpleTileTest {
         }
     }
 
+    @Test
+    public void testInterpolation() throws RuggedException {
+        SimpleTile tile = new SimpleTileFactory().createTile();
+        tile.setGeometry(0.0, 0.0, 1.0, 1.0, 50, 50);
+        tile.setElevation(20, 14,  91.0);
+        tile.setElevation(20, 15, 210.0);
+        tile.setElevation(21, 14, 162.0);
+        tile.setElevation(21, 15,  95.0);
+        Assert.assertEquals(150.5, tile.interpolateElevation(20.0, 14.5), 1.0e-10);
+        Assert.assertEquals(128.5, tile.interpolateElevation(21.0, 14.5), 1.0e-10);
+        Assert.assertEquals(146.1, tile.interpolateElevation(20.2, 14.5), 1.0e-10);
+    }
+
+    @Test
+    public void testInterpolationWithinTolerance() throws RuggedException {
+        SimpleTile tile = new SimpleTileFactory().createTile();
+        tile.setGeometry(0.0, 0.0, 1.0, 1.0, 2, 2);
+        tile.setElevation(0, 0,  91.0);
+        tile.setElevation(0, 1, 210.0);
+        tile.setElevation(1, 0, 162.0);
+        tile.setElevation(1, 1,  95.0);
+        // the following points are 1/16 pixel out of tile
+        Assert.assertEquals(151.875, tile.interpolateElevation(-0.0625,  0.5),    1.0e-10);
+        Assert.assertEquals(127.125, tile.interpolateElevation( 1.0625,  0.5),    1.0e-10);
+        Assert.assertEquals(124.875, tile.interpolateElevation( 0.5,    -0.0625), 1.0e-10);
+        Assert.assertEquals(154.125, tile.interpolateElevation( 0.5,     1.0625), 1.0e-10);
+    }
+
+    @Test
+    public void testInterpolationOutOfTolerance() throws RuggedException {
+        SimpleTile tile = new SimpleTileFactory().createTile();
+        tile.setGeometry(0.0, 0.0, 1.0, 1.0, 2, 2);
+        tile.setElevation(0, 0,  91.0);
+        tile.setElevation(0, 1, 210.0);
+        tile.setElevation(1, 0, 162.0);
+        tile.setElevation(1, 1,  95.0);
+        // the following points are 3/16 pixel out of tile
+        checkOutOfBound(-0.1875,  0.5,    tile);
+        checkOutOfBound( 1.1875,  0.5,    tile);
+        checkOutOfBound( 0.5,    -0.1875, tile);
+        checkOutOfBound( 0.5,     1.1875, tile);
+    }
+
+    private void checkOutOfBound(double latitude, double longitude, Tile tile) {
+        try {
+            tile.interpolateElevation(latitude, longitude);
+        } catch (RuggedException re) {
+            Assert.assertEquals(RuggedMessages.OUT_OF_TILE_ANGLES, re.getSpecifier());
+            Assert.assertEquals(FastMath.toDegrees(latitude),
+                                ((Double) re.getParts()[0]).doubleValue(),
+                                1.0e-10);
+            Assert.assertEquals(FastMath.toDegrees(longitude),
+                                ((Double) re.getParts()[1]).doubleValue(),
+                                1.0e-10);
+            Assert.assertEquals(FastMath.toDegrees(tile.getMinimumLatitude()),
+                                ((Double) re.getParts()[2]).doubleValue(),
+                                1.0e-10);
+            Assert.assertEquals(FastMath.toDegrees(tile.getMaximumLatitude()),
+                                ((Double) re.getParts()[3]).doubleValue(),
+                                1.0e-10);
+            Assert.assertEquals(FastMath.toDegrees(tile.getMinimumLongitude()),
+                                ((Double) re.getParts()[4]).doubleValue(),
+                                1.0e-10);
+            Assert.assertEquals(FastMath.toDegrees(tile.getMaximumLongitude()),
+                                ((Double) re.getParts()[5]).doubleValue(),
+                                1.0e-10);
+        }
+    }
+
 }
-- 
GitLab