From 0f831bc3a212c09f01ff24b4be1f92f5ce0cc70f Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Tue, 25 Mar 2014 20:02:57 +0100
Subject: [PATCH] Use point and direction rather than two points for pixel
intersection.
---
.../rugged/core/BasicScanAlgorithm.java | 2 +-
.../core/duvenhage/DuvenhageAlgorithm.java | 11 ++---
.../orekit/rugged/core/raster/SimpleTile.java | 48 +++++++++++--------
.../org/orekit/rugged/core/raster/Tile.java | 8 ++--
.../rugged/core/raster/SimpleTileTest.java | 33 ++++++++-----
5 files changed, 59 insertions(+), 43 deletions(-)
diff --git a/rugged-core/src/main/java/org/orekit/rugged/core/BasicScanAlgorithm.java b/rugged-core/src/main/java/org/orekit/rugged/core/BasicScanAlgorithm.java
index 49351cbd..8836cbe0 100644
--- a/rugged-core/src/main/java/org/orekit/rugged/core/BasicScanAlgorithm.java
+++ b/rugged-core/src/main/java/org/orekit/rugged/core/BasicScanAlgorithm.java
@@ -118,7 +118,7 @@ public class BasicScanAlgorithm implements IntersectionAlgorithm {
for (final SimpleTile tile : scannedTiles) {
for (int i = latitudeIndex(tile, minLatitude); i <= latitudeIndex(tile, maxLatitude); ++i) {
for (int j = longitudeIndex(tile, minLongitude); j <= longitudeIndex(tile, maxLongitude); ++j) {
- GeodeticPoint gp = tile.pixelIntersection(entryPoint, exitPoint, i, j);
+ GeodeticPoint gp = tile.pixelIntersection(entryPoint, ellipsoid.convertLos(entryPoint, los), i, j);
if (gp != null) {
final Vector3D point = ellipsoid.transform(gp);
final double dot = Vector3D.dotProduct(point.subtract(position), los);
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 b867f212..2f8aea9a 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
@@ -157,16 +157,13 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
final GeodeticPoint exit, final int exitLat, final int exitLon)
throws RuggedException, OrekitException {
- 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);
+ return tile.pixelIntersection(entry, ellipsoid.convertLos(entry, los), 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
@@ -320,9 +317,9 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
* @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;
- }
+ 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 {
diff --git a/rugged-core/src/main/java/org/orekit/rugged/core/raster/SimpleTile.java b/rugged-core/src/main/java/org/orekit/rugged/core/raster/SimpleTile.java
index cc69f80c..cae84330 100644
--- a/rugged-core/src/main/java/org/orekit/rugged/core/raster/SimpleTile.java
+++ b/rugged-core/src/main/java/org/orekit/rugged/core/raster/SimpleTile.java
@@ -16,6 +16,7 @@
*/
package org.orekit.rugged.core.raster;
+import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.Precision;
import org.orekit.bodies.GeodeticPoint;
@@ -244,7 +245,7 @@ public class SimpleTile implements Tile {
/** {@inheritDoc} */
@Override
- public GeodeticPoint pixelIntersection(final GeodeticPoint pA, final GeodeticPoint pB,
+ public GeodeticPoint pixelIntersection(final GeodeticPoint p, final Vector3D los,
final int latitudeIndex, final int longitudeIndex)
throws RuggedException {
@@ -257,15 +258,15 @@ public class SimpleTile implements Tile {
final double z11 = getElevationAtIndices(latitudeIndex + 1, longitudeIndex + 1);
// line-of-sight coordinates at close points
- final double dxA = (pA.getLongitude() - x00) / longitudeStep;
- final double dyA = (pA.getLatitude() - y00) / latitudeStep;
- final double dzA = pA.getAltitude();
- final double dxB = (pB.getLongitude() - x00) / longitudeStep;
- final double dyB = (pB.getLatitude() - y00) / latitudeStep;
- final double dzB = pB.getAltitude();
+ final double dxA = (p.getLongitude() - x00) / longitudeStep;
+ final double dyA = (p.getLatitude() - y00) / latitudeStep;
+ final double dzA = p.getAltitude();
+ final double dxB = dxA + los.getX() / longitudeStep;
+ final double dyB = dyA + los.getY() / latitudeStep;
+ final double dzB = dzA + los.getZ();
- // points along line-of-sight can be defined as a linear combination
- // between pA and pB depending on free variable t: p(t) = pA * (1 - t) + pB * t.
+ // points along line-of-sight can be defined as a linear progression
+ // along the line depending on free variable t: p(t) = p + t * los.
// As the point latitude and longitude are linear with respect to t,
// and as Digital Elevation Model is quadratic with respect to latitude
// and longitude, the altitude of DEM at same horizontal position as
@@ -307,8 +308,8 @@ public class SimpleTile implements Tile {
}
- final GeodeticPoint p1 = interpolate(t1, pA, dxA, dyA, pB, dxB, dyB);
- final GeodeticPoint p2 = interpolate(t2, pA, dxA, dyA, pB, dxB, dyB);
+ final GeodeticPoint p1 = interpolate(t1, p, dxA, dyA, los);
+ final GeodeticPoint p2 = interpolate(t2, p, dxA, dyA, los);
// select the first point along line-of-sight
if (p1 == null) {
@@ -321,23 +322,28 @@ public class SimpleTile implements Tile {
}
- /** Interpolate point on line.
- *
+ /** Interpolate point along a line.
+ * @param t abscissa along the line
+ * @param p start point
+ * @param dxP relative coordinate of the start point with respect to current pixel
+ * @param dyP relative coordinate of the start point with respect to current pixel
+ * @param los direction of the line-of-sight, in geodetic space
+ * @return interpolated point along the line
*/
- private GeodeticPoint interpolate(final double t,
- final GeodeticPoint pA, final double dxA, final double dyA,
- final GeodeticPoint pB, final double dxB, final double dyB) {
+ private GeodeticPoint interpolate(final double t, final GeodeticPoint p,
+ final double dxP, final double dyP,
+ final Vector3D los) {
if (Double.isInfinite(t)) {
return null;
}
- final double dx = dxA * (1 - t) + dxB * t;
- final double dy = dyA * (1 - t) + dyB * t;
+ final double dx = dxP + t * los.getX() / longitudeStep;
+ final double dy = dyP + t * los.getY() / latitudeStep;
if (dx >= 0 && dx <= 1 && dy >= 0 && dy <= 1) {
- return new GeodeticPoint(pA.getLatitude() * (1 - t) + pB.getLatitude() * t,
- pA.getLongitude() * (1 - t) + pB.getLongitude() * t,
- pA.getAltitude() * (1 - t) + pB.getAltitude() * t);
+ return new GeodeticPoint(p.getLatitude() + t * los.getY(),
+ p.getLongitude() + t * los.getX(),
+ p.getAltitude() + t * los.getZ());
} else {
return null;
}
diff --git a/rugged-core/src/main/java/org/orekit/rugged/core/raster/Tile.java b/rugged-core/src/main/java/org/orekit/rugged/core/raster/Tile.java
index 6a254e15..72ae22c3 100644
--- a/rugged-core/src/main/java/org/orekit/rugged/core/raster/Tile.java
+++ b/rugged-core/src/main/java/org/orekit/rugged/core/raster/Tile.java
@@ -16,6 +16,7 @@
*/
package org.orekit.rugged.core.raster;
+import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.api.RuggedException;
import org.orekit.rugged.api.UpdatableTile;
@@ -145,8 +146,9 @@ public interface Tile extends UpdatableTile {
throws RuggedException;
/** Find the intersection of a line-of-sight and a Digital Elevation Model pixel.
- * @param pA first point on the line (closer to satellite)
- * @param pB second point on the line (closer to ground)
+ * @param p point on the line
+ * @param los line-of-sight, in the topocentric frame (East, North, Zenith) of the point,
+ * scaled to match radians in the horizontal plane and meters along the vertical axis
* @param latitudeIndex latitude index of the Digital Elevation Model pixel
* @param longitudeIndex longitude index of the Digital Elevation Model pixel
* @return point corresponding to line-of-sight crossing the Digital Elevation Model surface
@@ -154,7 +156,7 @@ public interface Tile extends UpdatableTile {
* @exception RuggedException if intersection point cannot be computed
* @exception OrekitException if intersection point cannot be converted to geodetic coordinates
*/
- GeodeticPoint pixelIntersection(GeodeticPoint pA, GeodeticPoint pB,
+ GeodeticPoint pixelIntersection(GeodeticPoint p, Vector3D los,
int latitudeIndex, int longitudeIndex)
throws RuggedException;
diff --git a/rugged-core/src/test/java/org/orekit/rugged/core/raster/SimpleTileTest.java b/rugged-core/src/test/java/org/orekit/rugged/core/raster/SimpleTileTest.java
index a2867eca..7b68448a 100644
--- a/rugged-core/src/test/java/org/orekit/rugged/core/raster/SimpleTileTest.java
+++ b/rugged-core/src/test/java/org/orekit/rugged/core/raster/SimpleTileTest.java
@@ -16,6 +16,7 @@
*/
package org.orekit.rugged.core.raster;
+import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.apache.commons.math3.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
@@ -187,10 +188,10 @@ public class SimpleTileTest {
GeodeticPoint gpB = new GeodeticPoint(tile.getLatitudeAtIndex(20) + 0.7 * tile.getLatitudeStep(),
tile.getLongitudeAtIndex(14) + 0.9 * tile.getLongitudeStep(),
10.0);
- GeodeticPoint gpIAB = tile.pixelIntersection(gpA, gpB, 20, 14);
+ GeodeticPoint gpIAB = tile.pixelIntersection(gpA, los(gpA, gpB), 20, 14);
checkInLine(gpA, gpB, gpIAB);
checkOnTile(tile, gpIAB);
- GeodeticPoint gpIBA = tile.pixelIntersection(gpB, gpA, 20, 14);
+ GeodeticPoint gpIBA = tile.pixelIntersection(gpB, los(gpB, gpA), 20, 14);
checkInLine(gpA, gpB, gpIBA);
checkOnTile(tile, gpIBA);
@@ -218,10 +219,10 @@ public class SimpleTileTest {
// the line from gpA to gpB should traverse the DEM twice within the tile
// we use the points in the two different orders to retrieve both solutions
- GeodeticPoint gpIAB = tile.pixelIntersection(gpA, gpB, 20, 14);
+ GeodeticPoint gpIAB = tile.pixelIntersection(gpA, los(gpA, gpB), 20, 14);
checkInLine(gpA, gpB, gpIAB);
checkOnTile(tile, gpIAB);
- GeodeticPoint gpIBA = tile.pixelIntersection(gpB, gpA, 20, 14);
+ GeodeticPoint gpIBA = tile.pixelIntersection(gpB, los(gpB, gpA), 20, 14);
checkInLine(gpA, gpB, gpIBA);
checkOnTile(tile, gpIBA);
@@ -247,7 +248,7 @@ public class SimpleTileTest {
tile.getLongitudeAtIndex(14) + 0.9 * tile.getLongitudeStep(),
190.0);
- Assert.assertNull(tile.pixelIntersection(gpA, gpB, 20, 14));
+ Assert.assertNull(tile.pixelIntersection(gpA, los(gpA, gpB), 20, 14));
}
@@ -267,10 +268,10 @@ public class SimpleTileTest {
tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
20.0);
- GeodeticPoint gpIAB = tile.pixelIntersection(gpA, gpB, 0, 0);
+ GeodeticPoint gpIAB = tile.pixelIntersection(gpA, los(gpA, gpB), 0, 0);
checkInLine(gpA, gpB, gpIAB);
checkOnTile(tile, gpIAB);
- GeodeticPoint gpIBA = tile.pixelIntersection(gpB, gpA, 0, 0);
+ GeodeticPoint gpIBA = tile.pixelIntersection(gpB, los(gpB, gpA), 0, 0);
checkInLine(gpA, gpB, gpIBA);
checkOnTile(tile, gpIBA);
@@ -296,7 +297,7 @@ public class SimpleTileTest {
tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
55.0);
- Assert.assertNull(tile.pixelIntersection(gpA, gpB, 0, 0));
+ Assert.assertNull(tile.pixelIntersection(gpA, los(gpA, gpB), 0, 0));
}
@@ -316,7 +317,7 @@ public class SimpleTileTest {
tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
50.0);
- Assert.assertNull(tile.pixelIntersection(gpA, gpB, 0, 0));
+ Assert.assertNull(tile.pixelIntersection(gpA, los(gpA, gpB), 0, 0));
}
@@ -336,10 +337,10 @@ public class SimpleTileTest {
tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
37.5);
- GeodeticPoint gpIAB = tile.pixelIntersection(gpA, gpB, 0, 0);
+ GeodeticPoint gpIAB = tile.pixelIntersection(gpA, los(gpA, gpB), 0, 0);
checkInLine(gpA, gpB, gpIAB);
checkOnTile(tile, gpIAB);
- GeodeticPoint gpIBA = tile.pixelIntersection(gpB, gpA, 0, 0);
+ GeodeticPoint gpIBA = tile.pixelIntersection(gpB, los(gpB, gpA), 0, 0);
checkInLine(gpA, gpB, gpIBA);
checkOnTile(tile, gpIBA);
@@ -352,6 +353,16 @@ public class SimpleTileTest {
}
+ private Vector3D los(GeodeticPoint gpA, GeodeticPoint gpB) {
+ // this is a crude conversion into geodetic space
+ // intended *only* for the purposes of these tests
+ // it considers the geodetic space *is* perfectly Cartesian
+ // in the East, North, Zenith frame
+ return new Vector3D(gpB.getLongitude() - gpA.getLongitude(),
+ gpB.getLatitude() - gpA.getLatitude(),
+ gpB.getAltitude() - gpA.getAltitude());
+ }
+
private void checkOutOfBound(double latitude, double longitude, Tile tile) {
try {
tile.interpolateElevation(latitude, longitude);
--
GitLab