From 013115ea1e301b74384ae3d77b57d9aaecde9294 Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Thu, 18 Jun 2015 11:01:39 +0200
Subject: [PATCH] Added a protection against infinite loops in Duvenhage
algorithm.
---
.../duvenhage/DuvenhageAlgorithm.java | 134 +++++++++++++++---
src/site/xdoc/changes.xml | 3 +
.../org/orekit/rugged/api/RuggedTest.java | 2 +-
.../orekit/rugged/errors/DumpManagerTest.java | 2 +-
.../rugged/errors/DumpReplayerTest.java | 25 +++-
.../duvenhage/DuvenhageAlgorithmTest.java | 2 +-
.../resources/replay/replay-direct-loc-03.txt | 29 ++++
7 files changed, 176 insertions(+), 21 deletions(-)
create mode 100644 src/test/resources/replay/replay-direct-loc-03.txt
diff --git a/src/main/java/org/orekit/rugged/intersection/duvenhage/DuvenhageAlgorithm.java b/src/main/java/org/orekit/rugged/intersection/duvenhage/DuvenhageAlgorithm.java
index 33c13a4a..24e756b8 100644
--- a/src/main/java/org/orekit/rugged/intersection/duvenhage/DuvenhageAlgorithm.java
+++ b/src/main/java/org/orekit/rugged/intersection/duvenhage/DuvenhageAlgorithm.java
@@ -121,9 +121,9 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
final int exitLon = FastMath.max(0,
FastMath.min(tile.getLongitudeColumns() - 1,
tile.getFloorLongitudeIndex(exit.getPoint().getLongitude())));
- final NormalizedGeodeticPoint intersection = recurseIntersection(0, ellipsoid, position, los, tile,
- current, entryLat, entryLon,
- exit.getPoint(), exitLat, exitLon);
+ NormalizedGeodeticPoint intersection = recurseIntersection(0, ellipsoid, position, los, tile,
+ current, entryLat, entryLon,
+ exit.getPoint(), exitLat, exitLon);
if (intersection != null) {
// we have found the intersection
@@ -144,9 +144,19 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
}
} else {
+
// this should never happen
// we should have left the loop with an intersection point
- throw RuggedException.createInternalError(null);
+ // try a fallback non-recursive search
+ intersection = noRecurseIntersection(ellipsoid, position, los, tile,
+ current, entryLat, entryLon,
+ exitLat, exitLon);
+ if (intersection != null) {
+ return intersection;
+ } else {
+ throw RuggedException.createInternalError(null);
+ }
+
}
@@ -233,9 +243,10 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
throw RuggedException.createInternalError(null);
}
- if (entryLat == exitLat && entryLon == exitLon) {
- // we have narrowed the search down to a single Digital Elevation Model cell
- return tile.cellIntersection(entry, ellipsoid.convertLos(entry, los), exitLat, exitLon);
+ if (searchDomainSize(entryLat, entryLon, exitLat, exitLon) < 4) {
+ // we have narrowed the search down to a few cells
+ return noRecurseIntersection(ellipsoid, position, los, tile,
+ entry, entryLat, entryLon, exitLat, exitLon);
}
// find the deepest level in the min/max kd-tree at which entry and exit share a sub-tile
@@ -297,10 +308,18 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
if (inRange(crossingLonBefore, entryLon, exitLon)) {
// look for intersection
- final NormalizedGeodeticPoint intersection =
- recurseIntersection(depth + 1, ellipsoid, position, los, tile,
- previousGP, previousLat, previousLon,
- crossingGP, crossingLat, crossingLonBefore);
+ final NormalizedGeodeticPoint intersection;
+ if (searchDomainSize(previousLat, previousLon, crossingLat, crossingLonBefore) <
+ searchDomainSize(entryLat, entryLon, exitLat, exitLon)) {
+ intersection = recurseIntersection(depth + 1, ellipsoid, position, los, tile,
+ previousGP, previousLat, previousLon,
+ crossingGP, crossingLat, crossingLonBefore);
+ } else {
+ // we failed to reduce domain size, probably due to numerical problems
+ intersection = noRecurseIntersection(ellipsoid, position, los, tile,
+ previousGP, previousLat, previousLon,
+ crossingLat, crossingLonBefore);
+ }
if (intersection != null) {
return intersection;
}
@@ -359,9 +378,17 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
if (inRange(crossingLatBefore, entryLat, exitLat)) {
// look for intersection
- final NormalizedGeodeticPoint intersection = recurseIntersection(depth + 1, ellipsoid, position, los, tile,
- previousGP, previousLat, previousLon,
- crossingGP, crossingLatBefore, crossingLon);
+ final NormalizedGeodeticPoint intersection;
+ if (searchDomainSize(previousLat, previousLon, crossingLatBefore, crossingLon) <
+ searchDomainSize(entryLat, entryLon, exitLat, exitLon)) {
+ intersection = recurseIntersection(depth + 1, ellipsoid, position, los, tile,
+ previousGP, previousLat, previousLon,
+ crossingGP, crossingLatBefore, crossingLon);
+ } else {
+ intersection = noRecurseIntersection(ellipsoid, position, los, tile,
+ previousGP, previousLat, previousLon,
+ crossingLatBefore, crossingLon);
+ }
if (intersection != null) {
return intersection;
}
@@ -379,15 +406,88 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
if (inRange(previousLat, entryLat, exitLat) && inRange(previousLon, entryLon, exitLon)) {
// last part of the segment, up to exit point
- return recurseIntersection(depth + 1, ellipsoid, position, los, tile,
- previousGP, previousLat, previousLon,
- exit, exitLat, exitLon);
+ if (searchDomainSize(previousLat, previousLon, exitLat, exitLon) <
+ searchDomainSize(entryLat, entryLon, exitLat, exitLon)) {
+ return recurseIntersection(depth + 1, ellipsoid, position, los, tile,
+ previousGP, previousLat, previousLon,
+ exit, exitLat, exitLon);
+ } else {
+ return noRecurseIntersection(ellipsoid, position, los, tile,
+ previousGP, previousLat, previousLon,
+ exitLat, exitLon);
+ }
} else {
return null;
}
}
+ /** Compute intersection of line with Digital Elevation Model in a sub-tile, without recursion.
+ * @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 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 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 NormalizedGeodeticPoint noRecurseIntersection(final ExtendedEllipsoid ellipsoid,
+ final Vector3D position, final Vector3D los,
+ final MinMaxTreeTile tile,
+ final NormalizedGeodeticPoint entry,
+ final int entryLat, final int entryLon,
+ final int exitLat, final int exitLon)
+ throws OrekitException, RuggedException {
+
+ NormalizedGeodeticPoint intersectionGP = null;
+ double intersectionDot = Double.POSITIVE_INFINITY;
+ for (int i = FastMath.min(entryLat, exitLat); i <= FastMath.max(entryLat, exitLat); ++i) {
+ for (int j = FastMath.min(entryLon, exitLon); j <= FastMath.max(entryLon, exitLon); ++j) {
+ final NormalizedGeodeticPoint gp = tile.cellIntersection(entry, ellipsoid.convertLos(entry, los), i, j);
+ if (gp != null) {
+
+ // improve the point, by projecting it back on the 3D line, fixing the small body curvature at cell level
+ final Vector3D delta = ellipsoid.transform(gp).subtract(position);
+ final double s = Vector3D.dotProduct(delta, los) / los.getNormSq();
+ final GeodeticPoint projected = ellipsoid.transform(new Vector3D(1, position, s, los),
+ ellipsoid.getBodyFrame(), null);
+ final NormalizedGeodeticPoint gpImproved = tile.cellIntersection(projected, ellipsoid.convertLos(projected, los), i, j);
+
+ if (gpImproved != null) {
+ final Vector3D point = ellipsoid.transform(gpImproved);
+ final double dot = Vector3D.dotProduct(point.subtract(position), los);
+ if (dot < intersectionDot) {
+ intersectionGP = gpImproved;
+ intersectionDot = dot;
+ }
+ }
+
+ }
+ }
+ }
+
+ return intersectionGP;
+
+ }
+
+ /** Compute the size of a search domain.
+ * @param entryLat index to use for interpolating entry point elevation
+ * @param entryLon index to use for interpolating entry point elevation
+ * @param exitLat index to use for interpolating exit point elevation
+ * @param exitLon index to use for interpolating exit point elevation
+ * @return size of the search domain
+ */
+ private int searchDomainSize(final int entryLat, final int entryLon,
+ final int exitLat, final int exitLon) {
+ return (FastMath.abs(entryLat - exitLat) + 1) * (FastMath.abs(entryLon - exitLon) + 1);
+ }
+
/** Check if an index is inside a range.
* @param i index to check
* @param a first bound of the range (may be either below or above b)
diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml
index 8f9d7f45..9d1d111f 100644
--- a/src/site/xdoc/changes.xml
+++ b/src/site/xdoc/changes.xml
@@ -22,6 +22,9 @@
<body>
<release version="1.0" date="TBD"
description="TBD">
+ <action dev="luc" type="add" >
+ Added a protection against infinite loops in Duvenhage algorithm.
+ </action>
<action dev="luc" type="add" >
Added a CONSTANT_ELEVATION_OVER_ELLIPSOID algorithm, similar in spirit
to the IGNORE_DEM_USE_ELLIPSOID, but with a user-specified elevation
diff --git a/src/test/java/org/orekit/rugged/api/RuggedTest.java b/src/test/java/org/orekit/rugged/api/RuggedTest.java
index 00252890..853b4088 100644
--- a/src/test/java/org/orekit/rugged/api/RuggedTest.java
+++ b/src/test/java/org/orekit/rugged/api/RuggedTest.java
@@ -656,7 +656,7 @@ public class RuggedTest {
public void testInverseLocation()
throws RuggedException, OrekitException, URISyntaxException {
checkInverseLocation(2000, false, false, 3.0e-7, 5.0e-6);
- checkInverseLocation(2000, false, true, 6.0e-6, 2.0e-7);
+ checkInverseLocation(2000, false, true, 1.0e-5, 2.0e-7);
checkInverseLocation(2000, true, false, 4.0e-7, 4.0e-7);
checkInverseLocation(2000, true, true, 2.0e-5, 3.0e-7);
}
diff --git a/src/test/java/org/orekit/rugged/errors/DumpManagerTest.java b/src/test/java/org/orekit/rugged/errors/DumpManagerTest.java
index 10a12f37..46cb3d98 100644
--- a/src/test/java/org/orekit/rugged/errors/DumpManagerTest.java
+++ b/src/test/java/org/orekit/rugged/errors/DumpManagerTest.java
@@ -119,7 +119,7 @@ public class DumpManagerTest {
Assert.assertEquals(1, countSpan);
Assert.assertEquals(1, countTransform);
Assert.assertEquals(2, countDEMTile);
- Assert.assertEquals(887, countDEMCell);
+ Assert.assertEquals(812, countDEMCell);
Assert.assertEquals(400, countDirectLoc);
Assert.assertEquals(400, countDirectLocResult);
Assert.assertEquals(1, countSensor);
diff --git a/src/test/java/org/orekit/rugged/errors/DumpReplayerTest.java b/src/test/java/org/orekit/rugged/errors/DumpReplayerTest.java
index 572a74f2..31afa97d 100644
--- a/src/test/java/org/orekit/rugged/errors/DumpReplayerTest.java
+++ b/src/test/java/org/orekit/rugged/errors/DumpReplayerTest.java
@@ -51,7 +51,7 @@ public class DumpReplayerTest {
GeodeticPoint replayedGP = (GeodeticPoint) result.getReplayed();
double distance = Vector3D.distance(rugged.getEllipsoid().transform(expectedGP),
rugged.getEllipsoid().transform(replayedGP));
- Assert.assertEquals(0.0, distance, 1.0e-8);
+ Assert.assertEquals(0.0, distance, 6.0e-8);
}
}
@@ -79,6 +79,29 @@ public class DumpReplayerTest {
}
+ @Test
+ public void testDirectLoc03() throws URISyntaxException, IOException, OrekitException, RuggedException {
+
+ String orekitPath = getClass().getClassLoader().getResource("orekit-data").toURI().getPath();
+ DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(new File(orekitPath)));
+
+ String dumpPath = getClass().getClassLoader().getResource("replay/replay-direct-loc-03.txt").toURI().getPath();
+ DumpReplayer replayer = new DumpReplayer();
+ replayer.parse(new File(dumpPath));
+ Rugged rugged = replayer.createRugged();
+ DumpReplayer.Result[] results = replayer.execute(rugged);
+
+ Assert.assertEquals(1, results.length);
+ for (final DumpReplayer.Result result : results) {
+ GeodeticPoint expectedGP = (GeodeticPoint) result.getExpected();
+ GeodeticPoint replayedGP = (GeodeticPoint) result.getReplayed();
+ double distance = Vector3D.distance(rugged.getEllipsoid().transform(expectedGP),
+ rugged.getEllipsoid().transform(replayedGP));
+ Assert.assertEquals(0.0, distance, 1.0e-8);
+ }
+
+ }
+
@Test
public void testInverseLoc01() throws URISyntaxException, IOException, OrekitException, RuggedException {
diff --git a/src/test/java/org/orekit/rugged/intersection/duvenhage/DuvenhageAlgorithmTest.java b/src/test/java/org/orekit/rugged/intersection/duvenhage/DuvenhageAlgorithmTest.java
index 7f1e551f..a0a3fb51 100644
--- a/src/test/java/org/orekit/rugged/intersection/duvenhage/DuvenhageAlgorithmTest.java
+++ b/src/test/java/org/orekit/rugged/intersection/duvenhage/DuvenhageAlgorithmTest.java
@@ -119,7 +119,7 @@ public class DuvenhageAlgorithmTest extends AbstractAlgorithmTest {
algorithm.intersection(earth,
new Vector3D(-3041185.154503948, 6486750.132281409, -32335.022880173332),
new Vector3D(0.5660218606298548 , -0.8233939240951769, 0.040517885584811814));
- Assert.assertEquals(1164.339, gp.getAltitude(), 1.0e-3);
+ Assert.assertEquals(1164.35, gp.getAltitude(), 0.02);
}
@Test
diff --git a/src/test/resources/replay/replay-direct-loc-03.txt b/src/test/resources/replay/replay-direct-loc-03.txt
new file mode 100644
index 00000000..1643cd1a
--- /dev/null
+++ b/src/test/resources/replay/replay-direct-loc-03.txt
@@ -0,0 +1,29 @@
+# Rugged library dump file, created on 2015-06-17T12:40:59Z
+# all units are SI units (m, m/s, rad ...)
+direct location: date 2009-12-10T19:49:15.69409696600003Z position 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 los -2.222227735527040e-02 -9.103563670459193e-02 9.955996601239552e-01 lightTime false aberration false
+span: minDate 2009-12-10T19:43:11.00000000000000Z maxDate 2009-12-10T19:54:00.00000000000000Z tStep 1.000000000000000e-01 tolerance 1.000000000000000e+01 inertialFrame EME2000
+transform: index 3647 body r -9.892006943727071e-01 6.548538144006332e-05 4.928252940714298e-04 1.465665005635631e-01 Ω -7.256918583080753e-08 -1.041477941363816e-09 -7.292111444123439e-05 ΩDot -1.633367752479824e-16 1.747255489539140e-16 1.600518399853031e-19 spacecraft p 1.915468890219927e-01 -8.990556518547237e-01 -7.156923547128990e+06 v -2.488758510805678e-02 8.481932866150421e-02 -4.020837052430704e+00 a -2.066722296571628e-01 -2.090489712840085e+00 -7.003752909213540e-03 r -1.763451310531075e-01 -4.315699021299694e-01 8.655337192606597e-01 1.830333170551214e-01 Ω -8.737537450216512e-04 5.311070290168704e-04 2.103617306756347e-04 ΩDot 1.145446212396803e-07 1.919398959490098e-07 6.635983002829190e-08
+algorithm: DUVENHAGE
+ellipsoid: ae 6.378137000000000e+06 f 3.352810664747481e-03 frame ITRF_CIO_CONV_2010_SIMPLE_EOP
+DEM tile: t0 latMin 1.047197551196598e+00 latStep 1.454441043328608e-05 latRows 1201 lonMin -2.164208272472969e+00 lonStep 2.908882086657216e-05 lonCols 601
+DEM cell: t0 latIndex 1200 lonIndex 454 elevation 1.437147018278546e+02
+DEM cell: t0 latIndex 928 lonIndex 17 elevation 1.615747402326250e+03
+DEM cell: t0 latIndex 749 lonIndex 584 elevation 3.826208316815853e+02
+DEM cell: t0 latIndex 749 lonIndex 585 elevation 3.856180070818848e+02
+DEM cell: t0 latIndex 749 lonIndex 586 elevation 3.886151824821843e+02
+DEM cell: t0 latIndex 750 lonIndex 584 elevation 3.836203882299900e+02
+DEM cell: t0 latIndex 750 lonIndex 585 elevation 3.856175611585988e+02
+DEM cell: t0 latIndex 750 lonIndex 586 elevation 3.886147340872076e+02
+DEM cell: t0 latIndex 751 lonIndex 588 elevation 3.936086266060671e+02
+DEM cell: t0 latIndex 752 lonIndex 588 elevation 3.936081732677089e+02
+DEM cell: t0 latIndex 751 lonIndex 589 elevation 3.956057970629851e+02
+DEM cell: t0 latIndex 752 lonIndex 589 elevation 3.966053412529362e+02
+DEM cell: t0 latIndex 736 lonIndex 600 elevation 4.295819170687835e+02
+DEM cell: t0 latIndex 737 lonIndex 600 elevation 4.285814340701368e+02
+DEM cell: t0 latIndex 736 lonIndex 601 elevation 5.754850304128064e+02
+DEM cell: t0 latIndex 737 lonIndex 601 elevation 5.764867282588906e+02
+DEM cell: t0 latIndex 748 lonIndex 592 elevation 4.015986981091023e+02
+DEM cell: t0 latIndex 749 lonIndex 592 elevation 4.015982348839813e+02
+DEM cell: t0 latIndex 748 lonIndex 593 elevation 4.045958759810924e+02
+DEM cell: t0 latIndex 749 lonIndex 593 elevation 4.045954102842808e+02
+direct location result: latitude 1.058104370487310e+00 longitude -2.147166209870915e+00 elevation 3.882040377058536e+02
--
GitLab