From 62185321f01c2e8b65d62c5ed110f038cb8dc3ca Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Tue, 25 Mar 2014 14:30:33 +0100
Subject: [PATCH] First working version of Duvenhage algorithm!
---
.../core/duvenhage/DuvenhageAlgorithm.java | 18 +++---
.../rugged/core/duvenhage/MinMaxTreeTile.java | 59 +++++++++++--------
.../duvenhage/DuvenhageAlgorithmTest.java | 6 +-
.../core/duvenhage/MinMaxTreeTileTest.java | 4 +-
4 files changed, 52 insertions(+), 35 deletions(-)
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 a8e24776..b867f212 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
@@ -17,7 +17,6 @@
package org.orekit.rugged.core.duvenhage;
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
-import org.apache.commons.math3.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.errors.OrekitException;
import org.orekit.rugged.api.RuggedException;
@@ -152,20 +151,22 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
* @exception RuggedException if intersection cannot be found
* @exception OrekitException if points cannot be converted to geodetic coordinates
*/
- private GeodeticPoint recurseIntersection(final ExtendedEllipsoid ellipsoid,
- final Vector3D position, final Vector3D los,
- final MinMaxTreeTile tile,
+ 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 {
- if (FastMath.abs(entryLat - exitLat) < 1 && FastMath.abs(entryLon - exitLon) < 1) {
+ System.out.println(entryLat + " " + entryLon + " / " + exitLat + " " + exitLon);
+ if (entryLat == exitLat && entryLon == exitLon) {
// 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);
+ System.out.println("level = " + level + ", min = " + tile.getMinElevation(exitLat, exitLon, level) +
+ ", max = " + tile.getMaxElevation(exitLat, exitLon, level));
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
@@ -180,7 +181,9 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
// intersecting the boundary between level 0 sub-tiles
if (tile.isColumnMerging(level + 1)) {
// recurse through longitude crossings
- for (final int crossingLon : tile.getCrossedBoundaryColumns(previousLon, exitLon, level + 1)) {
+
+ int[] crossings = tile.getCrossedBoundaryColumns(previousLon, exitLon, level + 1);
+ for (final int crossingLon : crossings) {
// compute segment endpoints
final Vector3D crossingP = ellipsoid.pointAtLongitude(position, los,
@@ -208,7 +211,8 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
}
} else {
// recurse through latitude crossings
- for (final int crossingLat : tile.getCrossedBoundaryRows(previousLat, exitLat, level + 1)) {
+ int[] crossings = tile.getCrossedBoundaryRows(previousLat, exitLat, level + 1);
+ for (final int crossingLat : crossings) {
// compute segment endpoints
final Vector3D crossingP = ellipsoid.pointAtLatitude(position, los,
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 0118b77b..f5306264 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
@@ -106,6 +106,9 @@ public class MinMaxTreeTile extends SimpleTile {
// compute min/max trees
if (start.length > 0) {
+ // TODO: we need to pre-process pixels so the min/max of *interpolation* is used
+ // now, we consider a single corner in each pixel, but all four corners
+ // contributes to the min/max at interpolation stage
applyRecursively(minTree, start.length - 1, nbRows, nbCols, new Min(), elevations, 0);
applyRecursively(maxTree, start.length - 1, nbRows, nbCols, new Max(), elevations, 0);
}
@@ -209,22 +212,22 @@ public class MinMaxTreeTile extends SimpleTile {
* of these boundaries.
* </p>
* @param row1 starting row
- * @param row2 ending row (excluded)
+ * @param row2 ending row
* @param level tree level
- * @return indices of rows crossed at sub-tiles boundaries, in crossing order
+ * @return indices of rows crossed at sub-tiles boundaries, in crossing order,
+ * the endpoints <em>are</em> included (i.e. if {@code row1} or {@code row2} are
+ * boundary rows, they will be in returned array)
*/
public int[] getCrossedBoundaryRows(final int row1, final int row2, final int level) {
// number of rows in each sub-tile
- final int rows = 1 << ((start.length - level) / 2);
+ final int rows = 1 << ((start.length - level) / 2);
- if (row1 <= row2) {
- return buildCrossings((row1 + rows - 1) - ((row1 + rows - 1) % rows),
- row2, rows);
- } else {
- return buildCrossings(row1 -(row1 % rows),
- row2, -rows);
- }
+ // build the crossings in ascending order
+ final int min = FastMath.min(row1, row2);
+ final int max = FastMath.max(row1, row2) + 1;
+ return buildCrossings((min + rows - 1) - ((min + rows - 1) % rows), max, rows,
+ row1 <= row2);
}
@@ -237,30 +240,31 @@ public class MinMaxTreeTile extends SimpleTile {
* @param column1 starting column
* @param column2 ending column (excluded)
* @param level tree level
- * @return indices of columns crossed at sub-tiles boundaries, in crossing order
+ * @return indices of columns crossed at sub-tiles boundaries, in crossing order,
+ * the endpoints <em>are</em> included (i.e. if {@code column1} or {@code column2} are
+ * boundary columns, they will be in returned array)
*/
public int[] getCrossedBoundaryColumns(final int column1, final int column2, final int level) {
// number of columns in each sub-tile
- final int columns = 1 << ((start.length + 1 - level) / 2);;
+ final int columns = 1 << ((start.length + 1 - level) / 2);
- if (column1 <= column2) {
- return buildCrossings((column1 + columns - 1) - ((column1 + columns - 1) % columns),
- column2, columns);
- } else {
- return buildCrossings(column1 - (column1 % columns),
- column2, -columns);
- }
+ // build the crossings in ascending order
+ final int min = FastMath.min(column1, column2);
+ final int max = FastMath.max(column1, column2) + 1;
+ return buildCrossings((min + columns - 1) - ((min + columns - 1) % columns), max, columns,
+ column1 <= column2);
}
/** Build crossings arrays.
* @param begin begin crossing index
* @param end end crossing index (excluded, if equal to begin, the array is empty)
- * @param step crossing step (may be negative)
+ * @param step crossing step
+ * @param ascending if true, the crossings must be in ascending order
* @return indices of rows or columns crossed at sub-tiles boundaries, in crossing order
*/
- private int[] buildCrossings(final int begin, final int end, final int step) {
+ private int[] buildCrossings(final int begin, final int end, final int step, final boolean ascending) {
// allocate array
final int n = FastMath.max(0, (end - begin + step + ((step > 0) ? -1 : +1)) / step);
@@ -268,9 +272,16 @@ public class MinMaxTreeTile extends SimpleTile {
// fill it up
int crossing = begin;
- for (int i = 0; i < crossings.length; ++i) {
- crossings[i] = crossing;
- crossing += step;
+ if (ascending) {
+ for (int i = 0; i < crossings.length; ++i) {
+ crossings[i] = crossing;
+ crossing += step;
+ }
+ } else {
+ for (int i = 0; i < crossings.length; ++i) {
+ crossings[crossings.length - 1 - i] = crossing;
+ crossing += step;
+ }
}
return crossings;
diff --git a/rugged-core/src/test/java/org/orekit/rugged/core/duvenhage/DuvenhageAlgorithmTest.java b/rugged-core/src/test/java/org/orekit/rugged/core/duvenhage/DuvenhageAlgorithmTest.java
index 31caa691..eb9b75cd 100644
--- a/rugged-core/src/test/java/org/orekit/rugged/core/duvenhage/DuvenhageAlgorithmTest.java
+++ b/rugged-core/src/test/java/org/orekit/rugged/core/duvenhage/DuvenhageAlgorithmTest.java
@@ -65,6 +65,8 @@ public class DuvenhageAlgorithmTest {
FastMath.toRadians(1.0), 1201);
// test point approximately 1.6km North-North-West and 800 meters below volcano summit
+ // note that this test point is EXACTLY at a pixel corner, and even at corners of
+ // middle level (12 and above) sub-tiles
double latitude = FastMath.toRadians(13.27);
double longitude = FastMath.toRadians(123.68);
MinMaxTreeTile tile = new MinMaxTreeTileFactory().createTile();
@@ -126,7 +128,7 @@ public class DuvenhageAlgorithmTest {
GeodeticPoint result = duvenhage.intersection(earth, position, los);
Assert.assertEquals(groundGP.getLatitude(), result.getLatitude(), 1.0e-10);
Assert.assertEquals(groundGP.getLongitude(), result.getLongitude(), 1.0e-10);
- Assert.assertEquals(groundGP.getAltitude(), result.getAltitude(), 1.0e-10);
+ Assert.assertEquals(groundGP.getAltitude(), result.getAltitude(), 1.0e-9);
}
@@ -216,7 +218,7 @@ public class DuvenhageAlgorithmTest {
GeodeticPoint result = duvenhage.intersection(earth, position, los);
Assert.assertEquals(groundGP.getLatitude(), result.getLatitude(), 1.0e-10);
Assert.assertEquals(groundGP.getLongitude(), result.getLongitude(), 1.0e-10);
- Assert.assertEquals(groundGP.getAltitude(), result.getAltitude(), 1.0e-10);
+ Assert.assertEquals(groundGP.getAltitude(), result.getAltitude(), 1.0e-9);
}
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 280fa98c..72d30d26 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
@@ -254,8 +254,8 @@ 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 kS = FastMath.min(k1, k2);
+ int kE = FastMath.max(k1, k2) + 1;
int count = 0;
for (int k = kS; k < kE; ++k) {
if (k % n == 0) {
--
GitLab