Skip to content
Snippets Groups Projects
Commit 7c27c05d authored by Luc Maisonobe's avatar Luc Maisonobe
Browse files

Don't exclude endpoints when un-merging tiles.

parent e1791614
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/DuvenhageAlgorithm.java
+ 116
− 109
View file @ 7c27c05d
......@@ -16,9 +16,6 @@
*/
package org.orekit.rugged.core.duvenhage;
import java.util.ArrayList;
import java.util.List;
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.apache.commons.math3.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
......@@ -94,89 +91,42 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
// loop along the path
while (true) {
int currentLatIndex = tile.getLatitudeIndex(current.getLatitude());
int currentLonIndex = tile.getLongitudeIndex(current.getLongitude());
// find where line-of-sight exit tile
final LimitPoint exit = findExit(tile, ellipsoid, position, los);
final List<GeodeticPoint> lineOfSightQueue = new ArrayList<GeodeticPoint>();
lineOfSightQueue.add(exit.getPoint());
while (!lineOfSightQueue.isEmpty()) {
final GeodeticPoint next = lineOfSightQueue.remove(lineOfSightQueue.size() - 1);
final int nextLatIndex = tile.getLatitudeIndex(next.getLatitude());
final int nextLonIndex = tile.getLongitudeIndex(next.getLongitude());
if (FastMath.abs(currentLatIndex - nextLatIndex) <= 1 &&
FastMath.abs(currentLonIndex - nextLonIndex) <= 1) {
// we have narrowed the search down to a single Digital Elevation Model pixel
final GeodeticPoint intersection =
tile.pixelIntersection(current, next, nextLatIndex, nextLonIndex);
if (intersection != null) {
return intersection;
} else {
// no intersection on this pixel, we can proceed to next part of the line-of-sight
current = next;
}
} else {
// find the deepest level in the min/max kd-tree at which entry and exit share a sub-tile
final int level = tile.getMergeLevel(currentLatIndex, currentLonIndex, nextLatIndex, nextLonIndex);
if (level < 0) {
// introduce all intermediate points corresponding to the line-of-sight
// intersecting the boundary between level 0 sub-tiles
lineOfSightQueue.addAll(crossingPoints(ellipsoid, position, los,
tile, 0,
currentLatIndex, currentLonIndex,
nextLatIndex, nextLonIndex));
} else {
if (next.getAltitude() >= tile.getMaxElevation(nextLatIndex, nextLonIndex, 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
current = next;
} else {
// the line-of-sight segment has at least some undecided parts which may
// intersect the Digital Elevation Model, we need to refine the
// search by using a finer-grained level in the min/max kd-tree
// push the point back into the queue
lineOfSightQueue.add(next);
// introduce all intermediate points corresponding to the line-of-sight
// intersecting the boundary between finer sub-tiles as we go deeper
// in the tree
lineOfSightQueue.addAll(crossingPoints(ellipsoid, position, los,
tile, level,
currentLatIndex, currentLonIndex,
nextLatIndex, nextLonIndex));
// the current point remains the same
}
}
final GeodeticPoint intersection =
recurseIntersection(ellipsoid, position, los, tile,
current,
tile.getLatitudeIndex(current.getLatitude()),
tile.getLongitudeIndex(current.getLongitude()),
exit.getPoint(),
tile.getLatitudeIndex(exit.getPoint().getLatitude()),
tile.getLongitudeIndex(exit.getPoint().getLongitude()));
if (intersection != null) {
// we have found the intersection
return intersection;
} else if (exit.atSide()) {
// no intersection on this tile, we can proceed to next part of the line-of-sight
// select next tile after current point
final Vector3D forward = new Vector3D(1.0, ellipsoid.transform(exit.getPoint()), STEP, los);
current = ellipsoid.transform(forward, ellipsoid.getBodyFrame(), null);
tile = cache.getTile(current.getLatitude(), current.getLongitude());
if (tile.interpolateElevation(current.getLatitude(), current.getLongitude()) <= current.getAltitude()) {
// extremely rare case! The line-of-sight traversed the Digital Elevation Model
// during the very short forward step we used to move to next tile
// we consider this point to be OK
return current;
}
}
if (!exit.atSide()) {
} else {
// this should never happen
// we should have left the loop with an intersection point
throw RuggedException.createInternalError(null);
}
// select next tile after current point
final Vector3D forward = new Vector3D(1.0, ellipsoid.transform(current), STEP, los);
current = ellipsoid.transform(forward, ellipsoid.getBodyFrame(), null);
tile = cache.getTile(current.getLatitude(), current.getLongitude());
if (tile.interpolateElevation(current.getLatitude(), current.getLongitude()) <= current.getAltitude()) {
// extremely rare case! The line-of-sight traversed the Digital Elevation Model
// during the very short forward step we used to move to next tile
// we consider this point to be OK
return current;
}
}
......@@ -186,53 +136,110 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
}
}
/** Find the crossing points between sub-tiles.
* <p>
* When we go deeper in the min/max kd-tree, we get closer to individual pixels,
* or un-merge the merged sub-tiles. This un-merging implies finding the boundaries
* between sub-tiles that are merged at level l and not merged at level l+1.
* These boundaries are iso-latitude if the merge is a rows merging and are
* iso-longitude if the merge is a columns merging.
* </p>
/** Compute intersection of line with Digital Elevation Model in a sub-tile.
* @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 level merged level
* @param currentLatitude latitude index of current point (closer to observer)
* @param currentLongitude longitude index of current point (closer to observer)
* @param nextLatitude latitude index of next point (closer to ellipsoid)
* @param nextLongitude longitude index of next point (closer to ellipsoid)
* @return points corresponding to line-of-sight sub-tiles crossings, in
* <em>reversed</em> line-of-sight order
* @exception RuggedException if intersection point cannot be computed
* @exception OrekitException if intersection point cannot be converted to geodetic coordinates
* @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 exit line-of-sight exit point from the sub-tile
* @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 List<GeodeticPoint> crossingPoints(final ExtendedEllipsoid ellipsoid, final Vector3D position, final Vector3D los,
final MinMaxTreeTile tile, final int level,
final int currentLatitude, final int currentLongitude,
final int nextLatitude, final int nextLongitude)
private GeodeticPoint recurseIntersection(final ExtendedEllipsoid ellipsoid,
final Vector3D position, final Vector3D los,
final MinMaxTreeTile tile,
final GeodeticPoint entry, final int entryLat, final int entryLon,
final GeodeticPoint exit, final int exitLat, final int exitLon)
throws RuggedException, OrekitException {
final List<GeodeticPoint> crossings = new ArrayList<GeodeticPoint>();
if (FastMath.abs(entryLat - exitLat) < 1 && FastMath.abs(entryLon - exitLon) < 1) {
// we have narrowed the search down to a single Digital Elevation Model pixel
return tile.pixelIntersection(entry, exit, 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);
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
return null;
}
GeodeticPoint previousGP = entry;
int previousLat = entryLat;
int previousLon = entryLon;
// introduce all intermediate points corresponding to the line-of-sight
// intersecting the boundary between level 0 sub-tiles
if (tile.isColumnMerging(level + 1)) {
// sub-tiles at current level come from column merging at deeper level
for (final int longitudeIndex : tile.getCrossedBoundaryColumns(currentLongitude, nextLongitude, level + 1)) {
final double crossingLongitude = tile.getLongitudeAtIndex(longitudeIndex);
final Vector3D crossingPoint = ellipsoid.pointAtLongitude(position, los, crossingLongitude);
crossings.add(0, ellipsoid.transform(crossingPoint, ellipsoid.getBodyFrame(), null));
// recurse through longitude crossings
for (final int crossingLon : tile.getCrossedBoundaryColumns(previousLon, exitLon, level + 1)) {
// compute segment endpoints
final Vector3D crossingP = ellipsoid.pointAtLongitude(position, los,
tile.getLongitudeAtIndex(crossingLon));
final GeodeticPoint crossingGP = ellipsoid.transform(crossingP, ellipsoid.getBodyFrame(), null);
final int crossingLat = tile.getLatitudeIndex(crossingGP.getLatitude());
// adjust indices as the crossing point is by definition between the sub-tiles
final int crossingLonBefore = crossingLon - (entryLon <= exitLon ? 1 : 0);
final int crossingLonAfter = crossingLon - (entryLon <= exitLon ? 0 : 1);
// look for intersection
final GeodeticPoint intersection = recurseIntersection(ellipsoid, position, los, tile,
previousGP, previousLat, previousLon,
crossingGP, crossingLat, crossingLonBefore);
if (intersection != null) {
return intersection;
}
// prepare next segment
previousGP = crossingGP;
previousLat = crossingLat;
previousLon = crossingLonAfter;
}
} else {
// sub-tiles at current level come from row merging at deeper level
for (final int latitudeIndex : tile.getCrossedBoundaryRows(currentLatitude, nextLatitude, level + 1)) {
final double crossingLatitude = tile.getLatitudeAtIndex(latitudeIndex);
final Vector3D crossingPoint = ellipsoid.pointAtLatitude(position, los, crossingLatitude);
crossings.add(0, ellipsoid.transform(crossingPoint, ellipsoid.getBodyFrame(), null));
// recurse through latitude crossings
for (final int crossingLat : tile.getCrossedBoundaryRows(previousLat, exitLat, level + 1)) {
// compute segment endpoints
final Vector3D crossingP = ellipsoid.pointAtLatitude(position, los,
tile.getLatitudeAtIndex(crossingLat));
final GeodeticPoint crossingGP = ellipsoid.transform(crossingP, ellipsoid.getBodyFrame(), null);
final int crossingLon = tile.getLongitudeIndex(crossingGP.getLongitude());
// adjust indices as the crossing point is by definition between the sub-tiles
final int crossingLatBefore = crossingLat - (entryLat <= exitLat ? 1 : 0);
final int crossingLatAfter = crossingLat - (entryLat <= exitLat ? 0 : 1);
// look for intersection
final GeodeticPoint intersection = recurseIntersection(ellipsoid, position, los, tile,
previousGP, previousLat, previousLon,
crossingGP, crossingLatBefore, crossingLon);
if (intersection != null) {
return intersection;
}
// prepare next segment
previousGP = crossingGP;
previousLat = crossingLatAfter;
previousLon = crossingLon;
}
}
return crossings;
// last part of the segment, up to exit point
return recurseIntersection(ellipsoid, position, los, tile,
previousGP, previousLat, previousLon,
exit, exitLat, exitLon);
}
......
rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTile.java
+ 11
− 13
View file @ 7c27c05d
......@@ -209,10 +209,9 @@ public class MinMaxTreeTile extends SimpleTile {
* of these boundaries.
* </p>
* @param row1 starting row
* @param row2 ending row
* @param row2 ending row (excluded)
* @param level tree level
* @return indices of rows crossed at sub-tiles boundaries, in crossing order,
* excluding the start and end rows themselves even if they are at a sub-tile boundary
* @return indices of rows crossed at sub-tiles boundaries, in crossing order
*/
public int[] getCrossedBoundaryRows(final int row1, final int row2, final int level) {
......@@ -220,11 +219,11 @@ public class MinMaxTreeTile extends SimpleTile {
final int rows = 1 << ((start.length - level) / 2);
if (row1 <= row2) {
final int nextMultiple = row1 + rows - (row1 % rows);
return buildCrossings(nextMultiple, row2, rows);
return buildCrossings((row1 + rows - 1) - ((row1 + rows - 1) % rows),
row2, rows);
} else {
final int previousMultiple = row1 - 1 - ((row1 - 1) % rows);
return buildCrossings(previousMultiple, row2, -rows);
return buildCrossings(row1 -(row1 % rows),
row2, -rows);
}
}
......@@ -236,10 +235,9 @@ public class MinMaxTreeTile extends SimpleTile {
* of these boundaries.
* </p>
* @param column1 starting column
* @param column2 ending column
* @param column2 ending column (excluded)
* @param level tree level
* @return indices of columns crossed at sub-tiles boundaries, in crossing order
* excluding the start and end columns themselves even if they are at a sub-tile boundary
*/
public int[] getCrossedBoundaryColumns(final int column1, final int column2, final int level) {
......@@ -247,11 +245,11 @@ public class MinMaxTreeTile extends SimpleTile {
final int columns = 1 << ((start.length + 1 - level) / 2);;
if (column1 <= column2) {
final int nextMultiple = column1 + columns - (column1 % columns);
return buildCrossings(nextMultiple, column2, columns);
return buildCrossings((column1 + columns - 1) - ((column1 + columns - 1) % columns),
column2, columns);
} else {
final int previousMultiple = column1 - 1 - ((column1 - 1) % columns);
return buildCrossings(previousMultiple, column2, -columns);
return buildCrossings(column1 - (column1 % columns),
column2, -columns);
}
}
......
rugged-core/src/test/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java
+ 4
− 5
View file @ 7c27c05d
......@@ -206,9 +206,6 @@ public class MinMaxTreeTileTest {
for (int j2 = 0; j2 < nbColumns; ++j2) {
int[] crossings = tile.getCrossedBoundaryColumns(j1, j2, level);
int[] ref = multiples(j1, j2, subTileCols);
if (ref.length != crossings.length) {
crossings = tile.getCrossedBoundaryColumns(j1, j2, level);
}
Assert.assertArrayEquals(ref, crossings);
}
}
......@@ -257,8 +254,10 @@ public class MinMaxTreeTileTest {
// for testing purposes
// intentionally dumb way of counting multiples of n
int kS = (k1 <= k2) ? k1 : (k2 + 1);
int kE = (k1 <= k2) ? k2 : (k1 + 1);
int count = 0;
for (int k = FastMath.min(k1, k2) + 1; k < FastMath.max(k1, k2); ++k) {
for (int k = kS; k < kE; ++k) {
if (k % n == 0) {
++count;
}
......@@ -266,7 +265,7 @@ public class MinMaxTreeTileTest {
int[] multiples = new int[count];
int index = 0;
for (int k = FastMath.min(k1, k2) + 1; k < FastMath.max(k1, k2); ++k) {
for (int k = kS; k < kE; ++k) {
if (k % n == 0) {
multiples[index++] = k;
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment