diff --git a/src/changes/changes.xml b/src/changes/changes.xml index 3f4731dc63f6716c6d89e06d1a14af7ce5fb3697..75345e453f9de7f6067b0d2825924bd562226f81 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.1 (and Hipparchus 2.0). </action> diff --git a/src/main/java/org/orekit/rugged/errors/DumpReplayer.java b/src/main/java/org/orekit/rugged/errors/DumpReplayer.java index 2308472f7cd28059a2435011ed4781e192c1582e..4a02b3c56057a0237c21e7b1c04358dcd1ea788f 100644 --- a/src/main/java/org/orekit/rugged/errors/DumpReplayer.java +++ b/src/main/java/org/orekit/rugged/errors/DumpReplayer.java @@ -415,21 +415,7 @@ public class DumpReplayer { return rugged; - } catch (IOException ioe) { - throw new RuggedException(ioe, LocalizedCoreFormats.SIMPLE_MESSAGE, ioe.getLocalizedMessage()); - } catch (SecurityException e) { - // this should never happen - throw new RuggedInternalError(e); - } catch (NoSuchMethodException e) { - // this should never happen - throw new RuggedInternalError(e); - } catch (IllegalArgumentException e) { - // this should never happen - throw new RuggedInternalError(e); - } catch (IllegalAccessException e) { - // this should never happen - throw new RuggedInternalError(e); - } catch (InvocationTargetException e) { + } catch (IOException | NoSuchMethodException | IllegalAccessException | InvocationTargetException e) { // this should never happen throw new RuggedInternalError(e); } diff --git a/src/main/java/org/orekit/rugged/raster/SimpleTile.java b/src/main/java/org/orekit/rugged/raster/SimpleTile.java index 07e8dade5c703ed1d61abd65f3cccaf403d976b8..053c9993d5097dcccfca8a2b88084e11aa6bc14a 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 ce66a12607a7b4fe809a0d791882eb7aa6da90e0..d89506776c9e63ad38ce784981617d1e7a1409a7 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. @@ -257,7 +256,19 @@ public interface Tile extends UpdatableTile { double interpolateElevation(double latitude, double longitude); /** Find the intersection of a line-of-sight and a Digital Elevation Model cell. - * @param p point on the line + * <p> + * Beware that for continuity reasons, the point argument in {@code cellIntersection} is normalized + * with respect to other points used by the caller. This implies that the longitude may be + * outside of the [-π ; +π] interval (or the [0 ; 2π] interval, depending on the DEM). In particular, + * when a Line Of Sight crosses the antimeridian at ±π longitude, the library may call the + * {@code cellIntersection} method with a point having a longitude of -π-ε to ensure this continuity. + * As tiles are stored with longitude clipped to a some DEM specific interval (either [-π ; +π] or [0 ; 2π]), + * implementations MUST take care to clip the input point back to the tile interval using + * {@link org.hipparchus.util.MathUtils#normalizeAngle(double, double) MathUtils.normalizeAngle(p.getLongitude(), + * someLongitudeWithinTheTile)}. The output point normalization should also be made consistent with + * the current tile. + * </p> + * @param p point on the line (beware its longitude is <em>not</em> normalized with respect to tile) * @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 cell @@ -265,8 +276,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 a5b219ed3328bcb73862adc3823432512042db4d..64fcb92d2448eb532fce6f16cce803c3a9875c4a 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); }