diff --git a/rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/DuvenhageAlgorithm.java b/rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/DuvenhageAlgorithm.java
index 8950353ded0016967fa6c5f5dde0d07c6feb3291..a8e24776a75a1eff04a640cdd74fc679382e3b3a 100644
--- a/rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/DuvenhageAlgorithm.java
+++ b/rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/DuvenhageAlgorithm.java
@@ -16,9 +16,6 @@
*/
package org.orekit.rugged.core.duvenhage;
-import java.util.ArrayList;
-import java.util.List;
-
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.apache.commons.math3.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
@@ -94,89 +91,42 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
// loop along the path
while (true) {
- int currentLatIndex = tile.getLatitudeIndex(current.getLatitude());
- int currentLonIndex = tile.getLongitudeIndex(current.getLongitude());
-
// find where line-of-sight exit tile
final LimitPoint exit = findExit(tile, ellipsoid, position, los);
- final List<GeodeticPoint> lineOfSightQueue = new ArrayList<GeodeticPoint>();
- lineOfSightQueue.add(exit.getPoint());
-
- while (!lineOfSightQueue.isEmpty()) {
-
- final GeodeticPoint next = lineOfSightQueue.remove(lineOfSightQueue.size() - 1);
- final int nextLatIndex = tile.getLatitudeIndex(next.getLatitude());
- final int nextLonIndex = tile.getLongitudeIndex(next.getLongitude());
- if (FastMath.abs(currentLatIndex - nextLatIndex) <= 1 &&
- FastMath.abs(currentLonIndex - nextLonIndex) <= 1) {
-
- // we have narrowed the search down to a single Digital Elevation Model pixel
- final GeodeticPoint intersection =
- tile.pixelIntersection(current, next, nextLatIndex, nextLonIndex);
- if (intersection != null) {
- return intersection;
- } else {
- // no intersection on this pixel, we can proceed to next part of the line-of-sight
- current = next;
- }
-
- } else {
-
- // find the deepest level in the min/max kd-tree at which entry and exit share a sub-tile
- final int level = tile.getMergeLevel(currentLatIndex, currentLonIndex, nextLatIndex, nextLonIndex);
- if (level < 0) {
- // introduce all intermediate points corresponding to the line-of-sight
- // intersecting the boundary between level 0 sub-tiles
- lineOfSightQueue.addAll(crossingPoints(ellipsoid, position, los,
- tile, 0,
- currentLatIndex, currentLonIndex,
- nextLatIndex, nextLonIndex));
- } else {
- if (next.getAltitude() >= tile.getMaxElevation(nextLatIndex, nextLonIndex, level)) {
- // the line-of-sight 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 {
- // the line-of-sight segment has at least some undecided parts which may
- // intersect the Digital Elevation Model, we need to refine the
- // search by using a finer-grained level in the min/max kd-tree
-
- // push the point back into the queue
- lineOfSightQueue.add(next);
-
- // introduce all intermediate points corresponding to the line-of-sight
- // intersecting the boundary between finer sub-tiles as we go deeper
- // in the tree
- lineOfSightQueue.addAll(crossingPoints(ellipsoid, position, los,
- tile, level,
- currentLatIndex, currentLonIndex,
- nextLatIndex, nextLonIndex));
-
- // the current point remains the same
-
- }
- }
+ final GeodeticPoint intersection =
+ recurseIntersection(ellipsoid, position, los, tile,
+ current,
+ tile.getLatitudeIndex(current.getLatitude()),
+ tile.getLongitudeIndex(current.getLongitude()),
+ exit.getPoint(),
+ tile.getLatitudeIndex(exit.getPoint().getLatitude()),
+ tile.getLongitudeIndex(exit.getPoint().getLongitude()));
+
+ if (intersection != null) {
+ // we have found the intersection
+ return intersection;
+ } else if (exit.atSide()) {
+ // no intersection on this tile, we can proceed to next part of the line-of-sight
+
+ // select next tile after current point
+ final Vector3D forward = new Vector3D(1.0, ellipsoid.transform(exit.getPoint()), 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;
}
- }
-
- if (!exit.atSide()) {
+ } else {
// 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;
- }
}
@@ -186,53 +136,110 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
}
}
- /** Find the crossing points between sub-tiles.
- * <p>
- * When we go deeper in the min/max kd-tree, we get closer to individual pixels,
- * or un-merge the merged sub-tiles. This un-merging implies finding the boundaries
- * between sub-tiles that are merged at level l and not merged at level l+1.
- * These boundaries are iso-latitude if the merge is a rows merging and are
- * iso-longitude if the merge is a columns merging.
- * </p>
+ /** Compute intersection of line with Digital Elevation Model in a sub-tile.
* @param ellipsoid reference ellipsoid
* @param position pixel position in ellipsoid frame
* @param los pixel line-of-sight in ellipsoid frame
* @param tile Digital Elevation Model tile
- * @param level merged level
- * @param currentLatitude latitude index of current point (closer to observer)
- * @param currentLongitude longitude index of current point (closer to observer)
- * @param nextLatitude latitude index of next point (closer to ellipsoid)
- * @param nextLongitude longitude index of next point (closer to ellipsoid)
- * @return points corresponding to line-of-sight sub-tiles crossings, in
- * <em>reversed</em> line-of-sight order
- * @exception RuggedException if intersection point cannot be computed
- * @exception OrekitException if intersection point cannot be converted to geodetic coordinates
+ * @param entry line-of-sight entry point in the sub-tile
+ * @param entryLat index to use for interpolating entry point elevation
+ * @param entryLon index to use for interpolating entry point elevation
+ * @param exit line-of-sight exit point from the sub-tile
+ * @param exitLat index to use for interpolating exit point elevation
+ * @param exitLon index to use for interpolating exit point elevation
+ * @return point at which the line first enters ground, or null if does not enter
+ * ground in the search sub-tile
+ * @exception RuggedException if intersection cannot be found
+ * @exception OrekitException if points cannot be converted to geodetic coordinates
*/
- private List<GeodeticPoint> crossingPoints(final ExtendedEllipsoid ellipsoid, final Vector3D position, final Vector3D los,
- final MinMaxTreeTile tile, final int level,
- final int currentLatitude, final int currentLongitude,
- final int nextLatitude, final int nextLongitude)
+ private GeodeticPoint recurseIntersection(final ExtendedEllipsoid ellipsoid,
+ final Vector3D position, final Vector3D los,
+ final MinMaxTreeTile tile,
+ final GeodeticPoint entry, final int entryLat, final int entryLon,
+ final GeodeticPoint exit, final int exitLat, final int exitLon)
throws RuggedException, OrekitException {
- final List<GeodeticPoint> crossings = new ArrayList<GeodeticPoint>();
+ if (FastMath.abs(entryLat - exitLat) < 1 && FastMath.abs(entryLon - exitLon) < 1) {
+ // we have narrowed the search down to a single Digital Elevation Model pixel
+ return tile.pixelIntersection(entry, exit, exitLat, exitLon);
+ }
+
+ // find the deepest level in the min/max kd-tree at which entry and exit share a sub-tile
+ final int level = tile.getMergeLevel(entryLat, entryLon, exitLat, exitLon);
+ if (level >= 0 && exit.getAltitude() >= tile.getMaxElevation(exitLat, exitLon, level)) {
+ // the line-of-sight segment is fully above Digital Elevation Model
+ // we can safely reject it and proceed to next part of the line-of-sight
+ return null;
+ }
+ GeodeticPoint previousGP = entry;
+ int previousLat = entryLat;
+ int previousLon = entryLon;
+
+ // introduce all intermediate points corresponding to the line-of-sight
+ // intersecting the boundary between level 0 sub-tiles
if (tile.isColumnMerging(level + 1)) {
- // sub-tiles at current level come from column merging at deeper level
- for (final int longitudeIndex : tile.getCrossedBoundaryColumns(currentLongitude, nextLongitude, level + 1)) {
- final double crossingLongitude = tile.getLongitudeAtIndex(longitudeIndex);
- final Vector3D crossingPoint = ellipsoid.pointAtLongitude(position, los, crossingLongitude);
- crossings.add(0, ellipsoid.transform(crossingPoint, ellipsoid.getBodyFrame(), null));
+ // recurse through longitude crossings
+ for (final int crossingLon : tile.getCrossedBoundaryColumns(previousLon, exitLon, level + 1)) {
+
+ // compute segment endpoints
+ final Vector3D crossingP = ellipsoid.pointAtLongitude(position, los,
+ tile.getLongitudeAtIndex(crossingLon));
+ final GeodeticPoint crossingGP = ellipsoid.transform(crossingP, ellipsoid.getBodyFrame(), null);
+ final int crossingLat = tile.getLatitudeIndex(crossingGP.getLatitude());
+
+ // adjust indices as the crossing point is by definition between the sub-tiles
+ final int crossingLonBefore = crossingLon - (entryLon <= exitLon ? 1 : 0);
+ final int crossingLonAfter = crossingLon - (entryLon <= exitLon ? 0 : 1);
+
+ // look for intersection
+ final GeodeticPoint intersection = recurseIntersection(ellipsoid, position, los, tile,
+ previousGP, previousLat, previousLon,
+ crossingGP, crossingLat, crossingLonBefore);
+ if (intersection != null) {
+ return intersection;
+ }
+
+ // prepare next segment
+ previousGP = crossingGP;
+ previousLat = crossingLat;
+ previousLon = crossingLonAfter;
+
}
} else {
- // sub-tiles at current level come from row merging at deeper level
- for (final int latitudeIndex : tile.getCrossedBoundaryRows(currentLatitude, nextLatitude, level + 1)) {
- final double crossingLatitude = tile.getLatitudeAtIndex(latitudeIndex);
- final Vector3D crossingPoint = ellipsoid.pointAtLatitude(position, los, crossingLatitude);
- crossings.add(0, ellipsoid.transform(crossingPoint, ellipsoid.getBodyFrame(), null));
+ // recurse through latitude crossings
+ for (final int crossingLat : tile.getCrossedBoundaryRows(previousLat, exitLat, level + 1)) {
+
+ // compute segment endpoints
+ final Vector3D crossingP = ellipsoid.pointAtLatitude(position, los,
+ tile.getLatitudeAtIndex(crossingLat));
+ final GeodeticPoint crossingGP = ellipsoid.transform(crossingP, ellipsoid.getBodyFrame(), null);
+ final int crossingLon = tile.getLongitudeIndex(crossingGP.getLongitude());
+
+ // adjust indices as the crossing point is by definition between the sub-tiles
+ final int crossingLatBefore = crossingLat - (entryLat <= exitLat ? 1 : 0);
+ final int crossingLatAfter = crossingLat - (entryLat <= exitLat ? 0 : 1);
+
+ // look for intersection
+ final GeodeticPoint intersection = recurseIntersection(ellipsoid, position, los, tile,
+ previousGP, previousLat, previousLon,
+ crossingGP, crossingLatBefore, crossingLon);
+ if (intersection != null) {
+ return intersection;
+ }
+
+ // prepare next segment
+ previousGP = crossingGP;
+ previousLat = crossingLatAfter;
+ previousLon = crossingLon;
+
}
}
- return crossings;
+ // last part of the segment, up to exit point
+ return recurseIntersection(ellipsoid, position, los, tile,
+ previousGP, previousLat, previousLon,
+ exit, exitLat, exitLon);
}
diff --git a/rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTile.java b/rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTile.java
index 34ce31eda153f59c657540161390e5e0bc017723..0118b77b90e9c1daf6e338704d8a98d5067f98a8 100644
--- a/rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTile.java
+++ b/rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTile.java
@@ -209,10 +209,9 @@ public class MinMaxTreeTile extends SimpleTile {
* of these boundaries.
* </p>
* @param row1 starting row
- * @param row2 ending row
+ * @param row2 ending row (excluded)
* @param level tree level
- * @return indices of rows crossed at sub-tiles boundaries, in crossing order,
- * excluding the start and end rows themselves even if they are at a sub-tile boundary
+ * @return indices of rows crossed at sub-tiles boundaries, in crossing order
*/
public int[] getCrossedBoundaryRows(final int row1, final int row2, final int level) {
@@ -220,11 +219,11 @@ public class MinMaxTreeTile extends SimpleTile {
final int rows = 1 << ((start.length - level) / 2);
if (row1 <= row2) {
- final int nextMultiple = row1 + rows - (row1 % rows);
- return buildCrossings(nextMultiple, row2, rows);
+ return buildCrossings((row1 + rows - 1) - ((row1 + rows - 1) % rows),
+ row2, rows);
} else {
- final int previousMultiple = row1 - 1 - ((row1 - 1) % rows);
- return buildCrossings(previousMultiple, row2, -rows);
+ return buildCrossings(row1 -(row1 % rows),
+ row2, -rows);
}
}
@@ -236,10 +235,9 @@ public class MinMaxTreeTile extends SimpleTile {
* of these boundaries.
* </p>
* @param column1 starting column
- * @param column2 ending column
+ * @param column2 ending column (excluded)
* @param level tree level
* @return indices of columns crossed at sub-tiles boundaries, in crossing order
- * excluding the start and end columns themselves even if they are at a sub-tile boundary
*/
public int[] getCrossedBoundaryColumns(final int column1, final int column2, final int level) {
@@ -247,11 +245,11 @@ public class MinMaxTreeTile extends SimpleTile {
final int columns = 1 << ((start.length + 1 - level) / 2);;
if (column1 <= column2) {
- final int nextMultiple = column1 + columns - (column1 % columns);
- return buildCrossings(nextMultiple, column2, columns);
+ return buildCrossings((column1 + columns - 1) - ((column1 + columns - 1) % columns),
+ column2, columns);
} else {
- final int previousMultiple = column1 - 1 - ((column1 - 1) % columns);
- return buildCrossings(previousMultiple, column2, -columns);
+ return buildCrossings(column1 - (column1 % columns),
+ column2, -columns);
}
}
diff --git a/rugged-core/src/test/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java b/rugged-core/src/test/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java
index fbc66df073eb6480ab49c16a5815050de80851f5..280fa98c5a8c9d48800e5749157dd7902c34a5b1 100644
--- a/rugged-core/src/test/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java
+++ b/rugged-core/src/test/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java
@@ -206,9 +206,6 @@ public class MinMaxTreeTileTest {
for (int j2 = 0; j2 < nbColumns; ++j2) {
int[] crossings = tile.getCrossedBoundaryColumns(j1, j2, level);
int[] ref = multiples(j1, j2, subTileCols);
- if (ref.length != crossings.length) {
- crossings = tile.getCrossedBoundaryColumns(j1, j2, level);
- }
Assert.assertArrayEquals(ref, crossings);
}
}
@@ -257,8 +254,10 @@ public class MinMaxTreeTileTest {
// for testing purposes
// intentionally dumb way of counting multiples of n
+ int kS = (k1 <= k2) ? k1 : (k2 + 1);
+ int kE = (k1 <= k2) ? k2 : (k1 + 1);
int count = 0;
- for (int k = FastMath.min(k1, k2) + 1; k < FastMath.max(k1, k2); ++k) {
+ for (int k = kS; k < kE; ++k) {
if (k % n == 0) {
++count;
}
@@ -266,7 +265,7 @@ public class MinMaxTreeTileTest {
int[] multiples = new int[count];
int index = 0;
- for (int k = FastMath.min(k1, k2) + 1; k < FastMath.max(k1, k2); ++k) {
+ for (int k = kS; k < kE; ++k) {
if (k % n == 0) {
multiples[index++] = k;
}