From bfa8d8a7af76b550586c1fb95489a98b76d1500c Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Sat, 11 Dec 2021 10:38:32 +0100
Subject: [PATCH] Fixed longitude normalization issues.
Fixes #388
---
src/changes/changes.xml | 3 +
.../org/orekit/rugged/raster/SimpleTile.java | 33 +++--
.../java/org/orekit/rugged/raster/Tile.java | 5 +-
.../orekit/rugged/raster/SimpleTileTest.java | 136 +++++++++++-------
4 files changed, 107 insertions(+), 70 deletions(-)
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index fbc9ab78..af17c43f 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -21,6 +21,9 @@
</properties>
<body>
<release version="3.0" date="TBD" description="TBD">
+ <action dev="luc" type="update" issue="388">
+ Fixed longitude normalization issue with tiles.
+ </action>
<action dev="luc" type="update" issue="387">
Updated dependencies to Orekit 11.0.1 (and Hipparchus 2.0).
</action>
diff --git a/src/main/java/org/orekit/rugged/raster/SimpleTile.java b/src/main/java/org/orekit/rugged/raster/SimpleTile.java
index 07e8dade..053c9993 100644
--- a/src/main/java/org/orekit/rugged/raster/SimpleTile.java
+++ b/src/main/java/org/orekit/rugged/raster/SimpleTile.java
@@ -16,12 +16,11 @@
*/
package org.orekit.rugged.raster;
+import java.util.Arrays;
+
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.Precision;
-import java.util.Arrays;
-
-import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.errors.DumpManager;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.errors.RuggedMessages;
@@ -302,7 +301,7 @@ public class SimpleTile implements Tile {
/** {@inheritDoc} */
@Override
- public NormalizedGeodeticPoint cellIntersection(final GeodeticPoint p, final Vector3D los,
+ public NormalizedGeodeticPoint cellIntersection(final NormalizedGeodeticPoint p, final Vector3D los,
final int latitudeIndex, final int longitudeIndex) {
// ensure neighboring cells to not fall out of tile
@@ -317,10 +316,16 @@ public class SimpleTile implements Tile {
final double z10 = getElevationAtIndices(iLat, jLong + 1);
final double z11 = getElevationAtIndices(iLat + 1, jLong + 1);
+ // normalize back to tile coordinates
+ final NormalizedGeodeticPoint tileP = new NormalizedGeodeticPoint(p.getLatitude(),
+ p.getLongitude(),
+ p.getAltitude(),
+ x00);
+
// line-of-sight coordinates at close points
- final double dxA = (p.getLongitude() - x00) / longitudeStep;
- final double dyA = (p.getLatitude() - y00) / latitudeStep;
- final double dzA = p.getAltitude();
+ final double dxA = (tileP.getLongitude() - x00) / longitudeStep;
+ final double dyA = (tileP.getLatitude() - y00) / latitudeStep;
+ final double dzA = tileP.getAltitude();
final double dxB = dxA + los.getX() / longitudeStep;
final double dyB = dyA + los.getY() / latitudeStep;
final double dzB = dzA + los.getZ();
@@ -368,8 +373,8 @@ public class SimpleTile implements Tile {
}
- final NormalizedGeodeticPoint p1 = interpolate(t1, p, dxA, dyA, los, x00);
- final NormalizedGeodeticPoint p2 = interpolate(t2, p, dxA, dyA, los, x00);
+ final NormalizedGeodeticPoint p1 = interpolate(t1, tileP, dxA, dyA, los, x00);
+ final NormalizedGeodeticPoint p2 = interpolate(t2, tileP, dxA, dyA, los, x00);
// select the first point along line-of-sight
if (p1 == null) {
@@ -384,7 +389,7 @@ public class SimpleTile implements Tile {
/** Interpolate point along a line.
* @param t abscissa along the line
- * @param p start point
+ * @param tileP start point, normalized to tile area
* @param dxP relative coordinate of the start point with respect to current cell
* @param dyP relative coordinate of the start point with respect to current cell
* @param los direction of the line-of-sight, in geodetic space
@@ -392,7 +397,7 @@ public class SimpleTile implements Tile {
* be normalized between lc-Ï€ and lc+Ï€
* @return interpolated point along the line
*/
- private NormalizedGeodeticPoint interpolate(final double t, final GeodeticPoint p,
+ private NormalizedGeodeticPoint interpolate(final double t, final NormalizedGeodeticPoint tileP,
final double dxP, final double dyP,
final Vector3D los, final double centralLongitude) {
@@ -403,9 +408,9 @@ public class SimpleTile implements Tile {
final double dx = dxP + t * los.getX() / longitudeStep;
final double dy = dyP + t * los.getY() / latitudeStep;
if (dx >= -TOLERANCE && dx <= 1 + TOLERANCE && dy >= -TOLERANCE && dy <= 1 + TOLERANCE) {
- return new NormalizedGeodeticPoint(p.getLatitude() + t * los.getY(),
- p.getLongitude() + t * los.getX(),
- p.getAltitude() + t * los.getZ(),
+ return new NormalizedGeodeticPoint(tileP.getLatitude() + t * los.getY(),
+ tileP.getLongitude() + t * los.getX(),
+ tileP.getAltitude() + t * los.getZ(),
centralLongitude);
} else {
return null;
diff --git a/src/main/java/org/orekit/rugged/raster/Tile.java b/src/main/java/org/orekit/rugged/raster/Tile.java
index ce66a126..ccdb66dc 100644
--- a/src/main/java/org/orekit/rugged/raster/Tile.java
+++ b/src/main/java/org/orekit/rugged/raster/Tile.java
@@ -17,7 +17,6 @@
package org.orekit.rugged.raster;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
-import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.utils.NormalizedGeodeticPoint;
/** Interface representing a raster tile.
@@ -265,8 +264,8 @@ public interface Tile extends UpdatableTile {
* @return point corresponding to line-of-sight crossing the Digital Elevation Model surface
* if it lies within the cell, null otherwise
*/
- NormalizedGeodeticPoint cellIntersection(GeodeticPoint p, Vector3D los,
- int latitudeIndex, int longitudeIndex);
+ NormalizedGeodeticPoint cellIntersection(NormalizedGeodeticPoint p, Vector3D los,
+ int latitudeIndex, int longitudeIndex);
/** Check if a tile covers a ground point.
* @param latitude ground point latitude
diff --git a/src/test/java/org/orekit/rugged/raster/SimpleTileTest.java b/src/test/java/org/orekit/rugged/raster/SimpleTileTest.java
index a5b219ed..64fcb92d 100644
--- a/src/test/java/org/orekit/rugged/raster/SimpleTileTest.java
+++ b/src/test/java/org/orekit/rugged/raster/SimpleTileTest.java
@@ -18,12 +18,14 @@ package org.orekit.rugged.raster;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
+import org.hipparchus.util.MathUtils;
import org.junit.Assert;
import org.junit.Test;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.errors.RuggedMessages;
import org.orekit.rugged.raster.Tile.Location;
+import org.orekit.rugged.utils.NormalizedGeodeticPoint;
public class SimpleTileTest {
@@ -216,16 +218,44 @@ public class SimpleTileTest {
tile.setElevation(21, 14, 162.0);
tile.setElevation(21, 15, 95.0);
tile.tileUpdateCompleted();
- GeodeticPoint gpA = new GeodeticPoint(tile.getLatitudeAtIndex(20) + 0.1 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(14) + 0.2 * tile.getLongitudeStep(),
- 300.0);
- GeodeticPoint gpB = new GeodeticPoint(tile.getLatitudeAtIndex(20) + 0.7 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(14) + 0.9 * tile.getLongitudeStep(),
- 10.0);
- GeodeticPoint gpIAB = tile.cellIntersection(gpA, los(gpA, gpB), 20, 14);
+ NormalizedGeodeticPoint gpA = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(20) + 0.1 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(14) + 0.2 * tile.getLongitudeStep(),
+ 300.0, tile.getLongitudeAtIndex(14));
+ NormalizedGeodeticPoint gpB = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(20) + 0.7 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(14) + 0.9 * tile.getLongitudeStep(),
+ 10.0, tile.getLongitudeAtIndex(14));
+ NormalizedGeodeticPoint gpIAB = tile.cellIntersection(gpA, los(gpA, gpB), 20, 14);
checkInLine(gpA, gpB, gpIAB);
checkOnTile(tile, gpIAB);
- GeodeticPoint gpIBA = tile.cellIntersection(gpB, los(gpB, gpA), 20, 14);
+ NormalizedGeodeticPoint gpIBA = tile.cellIntersection(gpB, los(gpB, gpA), 20, 14);
+ checkInLine(gpA, gpB, gpIBA);
+ checkOnTile(tile, gpIBA);
+
+ Assert.assertEquals(gpIAB.getLatitude(), gpIBA.getLatitude(), 1.0e-10);
+ Assert.assertEquals(gpIAB.getLongitude(), gpIBA.getLongitude(), 1.0e-10);
+ Assert.assertEquals(gpIAB.getAltitude(), gpIBA.getAltitude(), 1.0e-10);
+
+ }
+
+ @Test
+ public void testCellIntersection2PiWrapping() {
+ SimpleTile tile = new SimpleTileFactory().createTile();
+ tile.setGeometry(0.0, 0.0, 0.025, 0.025, 50, 50);
+ tile.setElevation(20, 14, 91.0);
+ tile.setElevation(20, 15, 210.0);
+ tile.setElevation(21, 14, 162.0);
+ tile.setElevation(21, 15, 95.0);
+ tile.tileUpdateCompleted();
+ NormalizedGeodeticPoint gpA = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(20) + 0.1 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(14) + 0.2 * tile.getLongitudeStep(),
+ 300.0, +4 * FastMath.PI);
+ NormalizedGeodeticPoint gpB = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(20) + 0.7 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(14) + 0.9 * tile.getLongitudeStep(),
+ 10.0, +4 * FastMath.PI);
+ NormalizedGeodeticPoint gpIAB = tile.cellIntersection(gpA, los(gpA, gpB), 20, 14);
+ checkInLine(gpA, gpB, gpIAB);
+ checkOnTile(tile, gpIAB);
+ NormalizedGeodeticPoint gpIBA = tile.cellIntersection(gpB, los(gpB, gpA), 20, 14);
checkInLine(gpA, gpB, gpIBA);
checkOnTile(tile, gpIBA);
@@ -244,19 +274,19 @@ public class SimpleTileTest {
tile.setElevation(21, 14, 162.0);
tile.setElevation(21, 15, 95.0);
tile.tileUpdateCompleted();
- GeodeticPoint gpA = new GeodeticPoint(tile.getLatitudeAtIndex(20) + 0.1 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(14) + 0.2 * tile.getLongitudeStep(),
- 120.0);
- GeodeticPoint gpB = new GeodeticPoint(tile.getLatitudeAtIndex(20) + 0.7 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(14) + 0.9 * tile.getLongitudeStep(),
- 130.0);
+ NormalizedGeodeticPoint gpA = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(20) + 0.1 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(14) + 0.2 * tile.getLongitudeStep(),
+ 120.0, tile.getLongitudeAtIndex(14));
+ NormalizedGeodeticPoint gpB = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(20) + 0.7 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(14) + 0.9 * tile.getLongitudeStep(),
+ 130.0, tile.getLongitudeAtIndex(14));
// 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.cellIntersection(gpA, los(gpA, gpB), 20, 14);
+ NormalizedGeodeticPoint gpIAB = tile.cellIntersection(gpA, los(gpA, gpB), 20, 14);
checkInLine(gpA, gpB, gpIAB);
checkOnTile(tile, gpIAB);
- GeodeticPoint gpIBA = tile.cellIntersection(gpB, los(gpB, gpA), 20, 14);
+ NormalizedGeodeticPoint gpIBA = tile.cellIntersection(gpB, los(gpB, gpA), 20, 14);
checkInLine(gpA, gpB, gpIBA);
checkOnTile(tile, gpIBA);
@@ -275,12 +305,12 @@ public class SimpleTileTest {
tile.setElevation(21, 14, 162.0);
tile.setElevation(21, 15, 95.0);
tile.tileUpdateCompleted();
- GeodeticPoint gpA = new GeodeticPoint(tile.getLatitudeAtIndex(20) + 0.1 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(14) + 0.2 * tile.getLongitudeStep(),
- 180.0);
- GeodeticPoint gpB = new GeodeticPoint(tile.getLatitudeAtIndex(20) + 0.7 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(14) + 0.9 * tile.getLongitudeStep(),
- 190.0);
+ NormalizedGeodeticPoint gpA = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(20) + 0.1 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(14) + 0.2 * tile.getLongitudeStep(),
+ 180.0, tile.getLongitudeAtIndex(14));
+ NormalizedGeodeticPoint gpB = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(20) + 0.7 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(14) + 0.9 * tile.getLongitudeStep(),
+ 190.0, tile.getLongitudeAtIndex(14));
Assert.assertNull(tile.cellIntersection(gpA, los(gpA, gpB), 20, 14));
@@ -295,17 +325,17 @@ public class SimpleTileTest {
tile.setElevation(1, 0, 40.0);
tile.setElevation(1, 1, 40.0);
tile.tileUpdateCompleted();
- GeodeticPoint gpA = new GeodeticPoint(tile.getLatitudeAtIndex(0) + 0.25 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
- 50.0);
- GeodeticPoint gpB = new GeodeticPoint(tile.getLatitudeAtIndex(0) + 0.75 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
- 20.0);
-
- GeodeticPoint gpIAB = tile.cellIntersection(gpA, los(gpA, gpB), 0, 0);
+ NormalizedGeodeticPoint gpA = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(0) + 0.25 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
+ 50.0, tile.getLongitudeAtIndex(0));
+ NormalizedGeodeticPoint gpB = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(0) + 0.75 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
+ 20.0, tile.getLongitudeAtIndex(0));
+
+ NormalizedGeodeticPoint gpIAB = tile.cellIntersection(gpA, los(gpA, gpB), 0, 0);
checkInLine(gpA, gpB, gpIAB);
checkOnTile(tile, gpIAB);
- GeodeticPoint gpIBA = tile.cellIntersection(gpB, los(gpB, gpA), 0, 0);
+ NormalizedGeodeticPoint gpIBA = tile.cellIntersection(gpB, los(gpB, gpA), 0, 0);
checkInLine(gpA, gpB, gpIBA);
checkOnTile(tile, gpIBA);
@@ -324,12 +354,12 @@ public class SimpleTileTest {
tile.setElevation(1, 0, 40.0);
tile.setElevation(1, 1, 40.0);
tile.tileUpdateCompleted();
- GeodeticPoint gpA = new GeodeticPoint(tile.getLatitudeAtIndex(0) + 0.25 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
- 45.0);
- GeodeticPoint gpB = new GeodeticPoint(tile.getLatitudeAtIndex(0) + 0.75 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
- 55.0);
+ NormalizedGeodeticPoint gpA = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(0) + 0.25 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
+ 45.0, tile.getLongitudeAtIndex(0));
+ NormalizedGeodeticPoint gpB = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(0) + 0.75 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
+ 55.0, tile.getLongitudeAtIndex(0));
Assert.assertNull(tile.cellIntersection(gpA, los(gpA, gpB), 0, 0));
@@ -344,12 +374,12 @@ public class SimpleTileTest {
tile.setElevation(1, 0, 40.0);
tile.setElevation(1, 1, 40.0);
tile.tileUpdateCompleted();
- GeodeticPoint gpA = new GeodeticPoint(tile.getLatitudeAtIndex(0) + 0.25 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
- 45.0);
- GeodeticPoint gpB = new GeodeticPoint(tile.getLatitudeAtIndex(0) + 0.75 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
- 50.0);
+ NormalizedGeodeticPoint gpA = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(0) + 0.25 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
+ 45.0, tile.getLongitudeAtIndex(0));
+ NormalizedGeodeticPoint gpB = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(0) + 0.75 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
+ 50.0, tile.getLongitudeAtIndex(0));
Assert.assertNull(tile.cellIntersection(gpA, los(gpA, gpB), 0, 0));
@@ -364,17 +394,17 @@ public class SimpleTileTest {
tile.setElevation(1, 0, 40.0);
tile.setElevation(1, 1, 40.0);
tile.tileUpdateCompleted();
- GeodeticPoint gpA = new GeodeticPoint(tile.getLatitudeAtIndex(0) + 0.25 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
- 32.5);
- GeodeticPoint gpB = new GeodeticPoint(tile.getLatitudeAtIndex(0) + 0.75 * tile.getLatitudeStep(),
- tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
- 37.5);
-
- GeodeticPoint gpIAB = tile.cellIntersection(gpA, los(gpA, gpB), 0, 0);
+ NormalizedGeodeticPoint gpA = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(0) + 0.25 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
+ 32.5, tile.getLongitudeAtIndex(0));
+ NormalizedGeodeticPoint gpB = new NormalizedGeodeticPoint(tile.getLatitudeAtIndex(0) + 0.75 * tile.getLatitudeStep(),
+ tile.getLongitudeAtIndex(0) + 0.50 * tile.getLongitudeStep(),
+ 37.5, tile.getLongitudeAtIndex(0));
+
+ NormalizedGeodeticPoint gpIAB = tile.cellIntersection(gpA, los(gpA, gpB), 0, 0);
checkInLine(gpA, gpB, gpIAB);
checkOnTile(tile, gpIAB);
- GeodeticPoint gpIBA = tile.cellIntersection(gpB, los(gpB, gpA), 0, 0);
+ NormalizedGeodeticPoint gpIBA = tile.cellIntersection(gpB, los(gpB, gpA), 0, 0);
checkInLine(gpA, gpB, gpIBA);
checkOnTile(tile, gpIBA);
@@ -432,7 +462,7 @@ public class SimpleTileTest {
1.0e-10);
Assert.assertEquals(gpI.getLongitude(),
- gpA.getLongitude() * (1 - t) + gpB.getLongitude() * t,
+ MathUtils.normalizeAngle(gpA.getLongitude() * (1 - t) + gpB.getLongitude() * t, gpI.getLongitude()),
1.0e-10);
}
--
GitLab