started implementation of the Duvenhage algorithm.

rugged-api/src/main/java/org/orekit/rugged/api/RuggedMessages.java

@@ -52,12 +52,14 @@ public enum RuggedMessages implements Localizable {
 
     INTERNAL_ERROR("internal error, contact maintenance at {0}"),
     OUT_OF_TILE_INDICES("no data at indices [{0}, {1}], tile only covers from [0, 0] to [{2}, {3}] (inclusive)"),
+    OUT_OF_TILE_ANGLES("no data at latitude {0} and longitude {1}, tile covers only latitudes {2} to {3} and longitudes {4} to {5}"),
     UNINITIALIZED_CONTEXT("general context has not been initialized"),
     EMPTY_TILE("tile is empty: {0} ⨉ {1}"),
     UNKNOWN_SENSOR("unknown sensor {0}"),
     LINE_OF_SIGHT_NEVER_CROSSES_LATITUDE("line-of-sight never crosses latitude {0}"),
     LINE_OF_SIGHT_NEVER_CROSSES_LONGITUDE("line-of-sight never crosses longitude {0}"),
-    LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE("line-of-sight never crosses altitude {0}");
+    LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE("line-of-sight never crosses altitude {0}"),
+    DEM_ENTRY_POINT_IS_BEHIND_SPACECRAFT("line-of-sight enters the Digital Elevation Model behind spacecraft!");
 
     // CHECKSTYLE: resume JavadocVariable check
 

rugged-api/src/main/resources/assets/org/orekit/rugged/RuggedMessages_en.utf8

@@ -4,6 +4,9 @@ INTERNAL_ERROR = internal error, contact maintenance at {0}
 # no data at indices [{0}, {1}], tile only covers from [0, 0] to [{2}, {3}] (inclusive)
 OUT_OF_TILE_INDICES = no data at indices [{0}, {1}], tile only covers from [0, 0] to [{2}, {3}] (inclusive)
 
+# no data at latitude {0} and longitude {1}, tile covers only latitudes {2} to {3} and longitudes {4} to {5}
+OUT_OF_TILE_ANGLES = no data at latitude {0} and longitude {1}, tile covers only latitudes {2} to {3} and longitudes {4} to {5}
+
 # general context has not been initialized
 UNINITIALIZED_CONTEXT = general context has not been initialized
 
@@ -21,3 +24,6 @@ LINE_OF_SIGHT_NEVER_CROSSES_LONGITUDE = line-of-sight never crosses longitude {0
 
 # line never crosses altitude {0}
 LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE = line-of-sight never crosses altitude {0}
+
+# line-of-sight enters the Digital Elevation Model behind spacecraft!
+DEM_ENTRY_POINT_IS_BEHIND_SPACECRAFT = line-of-sight enters the Digital Elevation Model behind spacecraft!

rugged-api/src/main/resources/assets/org/orekit/rugged/RuggedMessages_fr.utf8

@@ -4,6 +4,9 @@ INTERNAL_ERROR = erreur interne, contactez la maintenance à {0}
 # no data at indices [{0}, {1}], tile only covers from [0, 0] to [{2}, {3}] (inclusive)
 OUT_OF_TILE_INDICES = aucune donnée aux indices [{0}, {1}], la tuile ne couvre que de [0, 0] à [{2}, {3}] inclus
 
+# no data at latitude {0} and longitude {1}, tile covers only latitudes {2} to {3} and longitudes {4} to {5}
+OUT_OF_TILE_ANGLES = aucune donnée à la latitude {0} et à la longitude {1}, la tuile ne couvre que les latitudes de {2} à {3} et les longitudes de {4} à {5}
+
 # general context has not been initialized
 UNINITIALIZED_CONTEXT = le contexte général n''a pas été initialisé
 
@@ -21,3 +24,6 @@ LINE_OF_SIGHT_NEVER_CROSSES_LONGITUDE = la ligne de visée ne franchit jamais la
 
 # line never crosses altitude {0}
 LINE_OF_SIGHT_NEVER_CROSSES_ALTITUDE = la ligne de visée ne franchit jamais l''altitude {0}
+
+# line-of-sight enters the Digital Elevation Model behind spacecraft!
+DEM_ENTRY_POINT_IS_BEHIND_SPACECRAFT = la ligne de visée entre dans le Modèle Numérique de Terrain derrière le satellite !

rugged-api/src/test/java/org/orekit/rugged/api/RuggedMessagesTest.java

@@ -29,7 +29,7 @@ public class RuggedMessagesTest {
 
     @Test
     public void testMessageNumber() {
-        Assert.assertEquals(8, RuggedMessages.values().length);
+        Assert.assertEquals(10, RuggedMessages.values().length);
     }
 
     @Test

rugged-core/src/main/java/org/orekit/rugged/core/dem/SimpleTile.java

@@ -183,16 +183,64 @@ public class SimpleTile implements Tile {
         return elevations[latitudeIndex * getLongitudeColumns() + longitudeIndex];
     }
 
+    /** {@inheritDoc} */
+    @Override
+    public double interpolateElevation(double latitude, double longitude)
+        throws RuggedException {
+
+        final double doubleLatitudeIndex  = getLatitudeIndex(latitude);
+        final double doubleLongitudeIndex = getLontitudeIndex(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)) {
+            throw new RuggedException(RuggedMessages.OUT_OF_TILE_ANGLES,
+                                      FastMath.toDegrees(latitude),
+                                      FastMath.toDegrees(longitude),
+                                      FastMath.toDegrees(getMinimumLatitude()),
+                                      FastMath.toDegrees(getMaximumLatitude()),
+                                      FastMath.toDegrees(getMinimumLongitude()),
+                                      FastMath.toDegrees(getMaximumLongitude()));
+        }
+
+        // bi-linear interpolation
+        final double dLat = doubleLatitudeIndex  - latitudeIndex;
+        final double dLon = doubleLongitudeIndex - longitudeIndex;
+        final double e00  = getElevationAtIndices(latitudeIndex,     longitudeIndex);
+        final double e01  = getElevationAtIndices(latitudeIndex,     longitudeIndex + 1);
+        final double e10  = getElevationAtIndices(latitudeIndex + 1, longitudeIndex);
+        final double e11  = getElevationAtIndices(latitudeIndex + 1, longitudeIndex + 1);
+        return (e00 * (1.0 - dLon) + dLon * e01) * (1.0 - dLat) +
+               (e10 * (1.0 - dLon) + dLon * e11) * dLat;
+
+    }
+
     /** {@inheritDoc} */
     @Override
     public int getLatitudeIndex(double latitude) {
-        return (int) FastMath.floor((latitude  - minLatitude)  / latitudeStep);
+        return (int) FastMath.floor(getDoubleLatitudeIndex(latitude));
     }
 
     /** {@inheritDoc} */
     @Override
     public int getLontitudeIndex(double longitude) {
-        return (int) FastMath.floor((longitude - minLongitude) / longitudeStep);
+        return (int) FastMath.floor(getDoubleLontitudeIndex(longitude));
+    }
+
+    /** Get the latitude index of a point.
+     * @param latitude geodetic latitude
+     * @return latirute index (it may lie outside of the tile!)
+     */
+    private double getDoubleLatitudeIndex(double latitude) {
+        return (latitude  - minLatitude)  / latitudeStep;
+    }
+
+    /** Get the longitude index of a point.
+     * @param longitude geodetic latitude
+     * @return longitude index (it may lie outside of the tile!)
+     */
+    private double getDoubleLontitudeIndex(double longitude) {
+        return (longitude - minLongitude) / longitudeStep;
     }
 
     /** {@inheritDoc} */

rugged-core/src/main/java/org/orekit/rugged/core/dem/Tile.java

@@ -114,6 +114,15 @@ public interface Tile extends UpdatableTile {
     double getElevationAtIndices(int latitudeIndex, int longitudeIndex)
         throws RuggedException;
 
+    /** Interpolate elevation.
+     * @param latitude ground point latitude
+     * @param longitude ground point longitude
+     * @return interpolated elevation
+     * @exception RuggedException if point does not lie within the tile
+     */
+    double interpolateElevation(double latitude, double longitude)
+        throws RuggedException;
+
     /** Check if a tile covers a ground point.
      * @param latitude ground point latitude
      * @param longitude ground point longitude

rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/DuvenhageAlgorithm.java

@@ -16,11 +16,18 @@
  */
 package org.orekit.rugged.core.duvenhage;
 
+import java.util.ArrayDeque;
+import java.util.Deque;
+
 import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
 import org.orekit.bodies.GeodeticPoint;
+import org.orekit.errors.OrekitException;
+import org.orekit.rugged.api.RuggedException;
+import org.orekit.rugged.api.RuggedMessages;
 import org.orekit.rugged.api.TileUpdater;
 import org.orekit.rugged.core.ExtendedEllipsoid;
 import org.orekit.rugged.core.dem.IntersectionAlgorithm;
+import org.orekit.rugged.core.dem.Tile;
 import org.orekit.rugged.core.dem.TilesCache;
 
 /** Digital Elevation Model intersection using Duvenhage's algorithm.
@@ -33,6 +40,9 @@ import org.orekit.rugged.core.dem.TilesCache;
  */
 public class DuvenhageAlgorithm implements IntersectionAlgorithm {
 
+    /** Step size when skipping from one tile to a neighbor one, in meters. */
+    private static final double STEP = 0.01;
+
     /** Cache for DEM tiles. */
     private TilesCache<MinMaxTreeTile> cache;
 
@@ -51,9 +61,218 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
     /** {@inheritDoc} */
     @Override
     public GeodeticPoint intersection(final ExtendedEllipsoid ellipsoid,
-                                      final Vector3D position, final Vector3D los) {
-        // TODO: compute intersection
-        return null;
+                                      final Vector3D position, final Vector3D los)
+        throws RuggedException {
+        try {
+
+            // compute intersection with ellipsoid
+            final Vector3D      p0  = ellipsoid.pointAtAltitude(position, los, 0.0);
+            final GeodeticPoint gp0 = ellipsoid.transform(p0, ellipsoid.getBodyFrame(), null);
+
+            // locate the entry tile along the line-of-sight
+            MinMaxTreeTile tile = cache.getTile(gp0.getLatitude(), gp0.getLongitude());
+
+            GeodeticPoint current = null;
+            while (current == null) {
+
+                // find where line-of-sight crosses tile max altitude
+                final Vector3D entryP = ellipsoid.pointAtAltitude(position, los, tile.getMaxElevation());
+                if (Vector3D.dotProduct(entryP.subtract(position), los) < 0) {
+                    // the entry point is behind spacecraft!
+                    throw new RuggedException(RuggedMessages.DEM_ENTRY_POINT_IS_BEHIND_SPACECRAFT);
+                }
+                current = ellipsoid.transform(entryP, ellipsoid.getBodyFrame(), null);
+
+                if (tile.getLocation(current.getLatitude(), current.getLongitude()) != Tile.Location.IN_TILE) {
+                    // the entry point is in another tile
+                    tile    = cache.getTile(current.getLatitude(), current.getLongitude());
+                    current = null;
+                } 
+
+            }
+
+            // loop along the path
+            while (true) {
+
+                int currentLatIndex = tile.getLatitudeIndex(current.getLatitude());
+                int currentLonIndex = tile.getLontitudeIndex(current.getLongitude());
+
+                // find where line-of-sight exit tile
+                final LimitPoint exit = findExit(tile, ellipsoid, position, los);
+                final Deque<GeodeticPoint> splitPointsQueue = new ArrayDeque<GeodeticPoint>();
+                splitPointsQueue.push(exit.getPoint());
+
+                while (!splitPointsQueue.isEmpty()) {
+
+                    final GeodeticPoint next = splitPointsQueue.pop();
+                    final int nextLatIndex   = tile.getLatitudeIndex(next.getLatitude());
+                    final int nextLonIndex   = tile.getLontitudeIndex(next.getLongitude());
+
+                    // find the largest level in the min/max kd-tree were entry and exit share a sub-tile
+                    int level = tile.getMergeLevel(currentLatIndex, currentLonIndex, nextLatIndex, nextLonIndex);
+                    if (level < 0) {
+                        // TODO: push intermediate points at sub-tiles boundaries on the queue
+                        throw RuggedException.createInternalError(null);
+                    }
+
+                    if (next.getAltitude() >= tile.getMaxElevation(nextLatIndex, nextLonIndex, level)) {
+                        // the line segment is fully above Digital Elevation Model
+                        // we can safely reject it and proceed to next part of the line-of-sight
+                        current = next;
+                    } else {
+                        // TODO: split line-of-sight
+                    }
+
+                }
+
+                if (!exit.atSide()) {
+                    // this should never happen
+                    // we should have left the loop with an intersection point
+                    throw RuggedException.createInternalError(null);                    
+                }
+
+                // select next tile after current point
+                final Vector3D forward = new Vector3D(1.0, ellipsoid.transform(current), STEP, los);
+                current = ellipsoid.transform(forward, ellipsoid.getBodyFrame(), null);
+                tile = cache.getTile(current.getLatitude(), current.getLongitude());
+                if (tile.interpolateElevation(current.getLatitude(), current.getLongitude()) <= current.getAltitude()) {
+                    // extremely rare case! The line-of-sight traversed the Digital Elevation Model
+                    // during the very short forward step we used to move to next tile
+                    // we consider this point to be OK
+                    return current;
+                }
+
+            }
+
+
+        } catch (OrekitException oe) {
+            throw new RuggedException(oe, oe.getSpecifier(), oe.getParts());
+        }
+    }
+
+    /** Compute a line-of-sight exit point from a tile.
+     * @param tile tile to consider
+     * @param ellipsoid reference ellipsoid
+     * @param position pixel position in ellipsoid frame
+     * @param los pixel line-of-sight in ellipsoid frame
+     * @return exit point
+     * @exception RuggedException if exit point cannot be found
+     * @exception OrekitException if geodetic coordinates cannot be computed
+     */
+    private LimitPoint findExit(final Tile tile, final ExtendedEllipsoid ellipsoid,
+                                final Vector3D position, final Vector3D los)
+        throws RuggedException, OrekitException {
+
+        // look for an exit at bottom
+        Vector3D exitP = ellipsoid.pointAtAltitude(position, los, tile.getMinElevation());
+        GeodeticPoint exitGP = ellipsoid.transform(exitP, ellipsoid.getBodyFrame(), null);
+
+        switch (tile.getLocation(exitGP.getLatitude(), exitGP.getLongitude())) {
+            case SOUTH_WEST :
+                return new LimitPoint(ellipsoid,
+                                      selectClosest(ellipsoid.pointAtLatitude(position,  los, tile.getMinimumLatitude()),
+                                                    ellipsoid.pointAtLongitude(position, los, tile.getMinimumLongitude()),
+                                                    position),
+                                      true);
+            case WEST :
+                return new LimitPoint(ellipsoid,
+                                      ellipsoid.pointAtLongitude(position, los, tile.getMinimumLongitude()),
+                                      true);
+            case NORTH_WEST:
+                return new LimitPoint(ellipsoid,
+                                      selectClosest(ellipsoid.pointAtLatitude(position,  los, tile.getMaximumLatitude()),
+                                                    ellipsoid.pointAtLongitude(position, los, tile.getMinimumLongitude()),
+                                                    position),
+                                      true);
+            case NORTH :
+                return new LimitPoint(ellipsoid,
+                                      ellipsoid.pointAtLatitude(position, los, tile.getMaximumLatitude()),
+                                      true);
+            case NORTH_EAST :
+                return new LimitPoint(ellipsoid,
+                                      selectClosest(ellipsoid.pointAtLatitude(position,  los, tile.getMaximumLatitude()),
+                                                    ellipsoid.pointAtLongitude(position, los, tile.getMaximumLongitude()),
+                                                    position),
+                                      true);
+            case EAST :
+                return new LimitPoint(ellipsoid,
+                                      ellipsoid.pointAtLongitude(position, los, tile.getMaximumLongitude()),
+                                      true);
+            case SOUTH_EAST :
+                return new LimitPoint(ellipsoid,
+                                      selectClosest(ellipsoid.pointAtLatitude(position,  los, tile.getMinimumLatitude()),
+                                                    ellipsoid.pointAtLongitude(position, los, tile.getMaximumLongitude()),
+                                                    position),
+                                      true);
+            case SOUTH :
+                return new LimitPoint(ellipsoid,
+                                      ellipsoid.pointAtLatitude(position, los, tile.getMinimumLatitude()),
+                                      true);
+            case IN_TILE :
+                return new LimitPoint(exitGP, false);
+            default :
+                // this should never happen
+                throw RuggedException.createInternalError(null);
+        }
+        
+    }
+
+    /** Select point closest to line-of-sight start.
+     * @param p1 first point to consider
+     * @param p2 second point to consider
+     * @param position pixel position in ellipsoid frame
+     * @return either p1 or p2, depending on which is closest to position
+     */
+   private Vector3D selectClosest(Vector3D p1, Vector3D p2, Vector3D position) {
+       return Vector3D.distance(p1, position) <= Vector3D.distance(p2, position) ? p1 : p2;
+   }
+
+    /** Point at tile boundary. */
+    private static class LimitPoint {
+
+        /** Coordinates. */
+        private final GeodeticPoint point;
+
+        /** Limit status. */
+        private final boolean side;
+
+        /** Simple constructor.
+         * @param cartesian point cartesian
+         * @param ellipsoid reference ellipsoid
+         * @param side if true, the point is on a side limit, otherwise
+         * it is on a top/bottom limit
+         * @exception OrekitException if geodetic coordinates cannot be computed
+         */
+        public LimitPoint(final ExtendedEllipsoid ellipsoid, final Vector3D cartesian, final boolean side)
+            throws OrekitException {
+            this(ellipsoid.transform(cartesian, ellipsoid.getBodyFrame(), null), side);
+        }
+
+        /** Simple constructor.
+         * @param point coordinates
+         * @param side if true, the point is on a side limit, otherwise
+         * it is on a top/bottom limit
+         */
+        public LimitPoint(final GeodeticPoint point, final boolean side) {
+            this.point = point;
+            this.side  = side;
+        }
+
+        /** Get the point coordinates.
+         * @return point coordinates
+         */
+        public GeodeticPoint getPoint() {
+            return point;
+        }
+
+        /** Check if point is on the side of a tile.
+         * @return true if the point is on a side limit, otherwise
+         * it is on a top/bottom limit
+         */
+        public boolean atSide() {
+            return side;
+        }
+
     }
 
 }