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 8369ee2cc8751b9f385c359de5e2aebd38523e6b..cc69f80c9023ceab9f00c2d6731c516975d675c2 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
@@ -17,6 +17,7 @@
package org.orekit.rugged.core.raster;
import org.apache.commons.math3.util.FastMath;
+import org.apache.commons.math3.util.Precision;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.api.RuggedException;
import org.orekit.rugged.api.RuggedMessages;
@@ -284,39 +285,31 @@ public class SimpleTile implements Tile {
final double b = v + dzA - dzB;
final double c = w - dzA;
- // solve the quadratic equation
- final double b2 = b * b;
- final double fac = 4 * a * c;
- if (b2 < fac) {
- // no intersection at all
- return null;
- }
- final double s = FastMath.sqrt(b2 - fac);
- final double t1 = (b < 0) ? (s - b) / (2 * a) : -2 * c / (b + s);
- final double t2 = c / (a * t1);
-
- final double dx1 = dxA * (1 - t1) + dxB * t1;
- final double dy1 = dyA * (1 - t1) + dyB * t1;
- final GeodeticPoint p1;
- if (dx1 >= 0 && dx1 <= 1 && dy1 >= 0 && dy1 <= 1) {
- p1 = new GeodeticPoint(pA.getLatitude() * (1 - t1) + pB.getLatitude() * t1,
- pA.getLongitude() * (1 - t1) + pB.getLongitude() * t1,
- pA.getAltitude() * (1 - t1) + pB.getAltitude() * t1);
+ // solve the equation
+ final double t1;
+ final double t2;
+ if (FastMath.abs(a) <= Precision.EPSILON * FastMath.abs(c)) {
+ // the equation degenerates to a linear (or constant) equation
+ double t = -c / b;
+ t1 = Double.isNaN(t) ? 0.0 : t;
+ t2 = Double.POSITIVE_INFINITY;
} else {
- p1 = null;
- }
+ // the equation is quadratic
+ final double b2 = b * b;
+ final double fac = 4 * a * c;
+ if (b2 < fac) {
+ // no intersection at all
+ return null;
+ }
+ final double s = FastMath.sqrt(b2 - fac);
+ t1 = (b < 0) ? (s - b) / (2 * a) : -2 * c / (b + s);
+ t2 = c / (a * t1);
- final double dx2 = dxA * (1 - t2) + dxB * t2;
- final double dy2 = dyA * (1 - t2) + dyB * t2;
- final GeodeticPoint p2;
- if (dx2 >= 0 && dx2 <= 1 && dy2 >= 0 && dy2 <= 1) {
- p2 = new GeodeticPoint(pA.getLatitude() * (1 - t2) + pB.getLatitude() * t2,
- pA.getLongitude() * (1 - t2) + pB.getLongitude() * t2,
- pA.getAltitude() * (1 - t2) + pB.getAltitude() * t2);
- } else {
- p2 = null;
}
+ final GeodeticPoint p1 = interpolate(t1, pA, dxA, dyA, pB, dxB, dyB);
+ final GeodeticPoint p2 = interpolate(t2, pA, dxA, dyA, pB, dxB, dyB);
+
// select the first point along line-of-sight
if (p1 == null) {
return p2;
@@ -328,6 +321,29 @@ public class SimpleTile implements Tile {
}
+ /** Interpolate point on 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) {
+
+ if (Double.isInfinite(t)) {
+ return null;
+ }
+
+ final double dx = dxA * (1 - t) + dxB * t;
+ final double dy = dyA * (1 - t) + dyB * t;
+ 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);
+ } else {
+ return null;
+ }
+
+ }
+
/** {@inheritDoc} */
@Override
public int getLatitudeIndex(double latitude) {
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 363d39d2e7b23b8758768158ccdd29172752f4bb..a2867eca00fb53134c9e69cee489c09b6d825973 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
@@ -251,6 +251,107 @@ public class SimpleTileTest {
}
+ @Test
+ public void testPixelIntersectionLinearOnly() throws RuggedException {
+ SimpleTile tile = new SimpleTileFactory().createTile();
+ tile.setGeometry(0.0, 0.0, 0.025, 0.025, 50, 50);
+ tile.setElevation(0, 0, 30.0);
+ tile.setElevation(0, 1, 30.0);
+ 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.pixelIntersection(gpA, gpB, 0, 0);
+ checkInLine(gpA, gpB, gpIAB);
+ checkOnTile(tile, gpIAB);
+ GeodeticPoint gpIBA = tile.pixelIntersection(gpB, gpA, 0, 0);
+ 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 testPixelIntersectionLinearIntersectionOutside() throws RuggedException {
+ SimpleTile tile = new SimpleTileFactory().createTile();
+ tile.setGeometry(0.0, 0.0, 0.025, 0.025, 50, 50);
+ tile.setElevation(0, 0, 30.0);
+ tile.setElevation(0, 1, 30.0);
+ 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);
+
+ Assert.assertNull(tile.pixelIntersection(gpA, gpB, 0, 0));
+
+ }
+
+ @Test
+ public void testPixelIntersectionLinearNoIntersection() throws RuggedException {
+ SimpleTile tile = new SimpleTileFactory().createTile();
+ tile.setGeometry(0.0, 0.0, 0.025, 0.025, 50, 50);
+ tile.setElevation(0, 0, 30.0);
+ tile.setElevation(0, 1, 30.0);
+ 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);
+
+ Assert.assertNull(tile.pixelIntersection(gpA, gpB, 0, 0));
+
+ }
+
+ @Test
+ public void testPixelIntersectionConstant0() throws RuggedException {
+ SimpleTile tile = new SimpleTileFactory().createTile();
+ tile.setGeometry(0.0, 0.0, 0.025, 0.025, 50, 50);
+ tile.setElevation(0, 0, 30.0);
+ tile.setElevation(0, 1, 30.0);
+ 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.pixelIntersection(gpA, gpB, 0, 0);
+ checkInLine(gpA, gpB, gpIAB);
+ checkOnTile(tile, gpIAB);
+ GeodeticPoint gpIBA = tile.pixelIntersection(gpB, gpA, 0, 0);
+ checkInLine(gpA, gpB, gpIBA);
+ checkOnTile(tile, gpIBA);
+
+ Assert.assertEquals(gpIAB.getLatitude(), gpA.getLatitude(), 1.0e-10);
+ Assert.assertEquals(gpIAB.getLongitude(), gpA.getLongitude(), 1.0e-10);
+ Assert.assertEquals(gpIAB.getAltitude(), gpA.getAltitude(), 1.0e-10);
+ Assert.assertEquals(gpIBA.getLatitude(), gpB.getLatitude(), 1.0e-10);
+ Assert.assertEquals(gpIBA.getLongitude(), gpB.getLongitude(), 1.0e-10);
+ Assert.assertEquals(gpIBA.getAltitude(), gpB.getAltitude(), 1.0e-10);
+
+ }
+
private void checkOutOfBound(double latitude, double longitude, Tile tile) {
try {
tile.interpolateElevation(latitude, longitude);