From 958739ab7868044e217b722ca75f6e2f55beac84 Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Tue, 18 Mar 2014 16:48:24 +0100
Subject: [PATCH] Added line-of-sight splitting at sub-tiles boundaries.
---
.../core/duvenhage/DuvenhageAlgorithm.java | 103 +++++++++----
.../rugged/core/duvenhage/MinMaxTreeTile.java | 141 ++++++++++++++----
.../core/duvenhage/MinMaxTreeTileTest.java | 59 ++++----
3 files changed, 209 insertions(+), 94 deletions(-)
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 b786117b..a6e805f2 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
@@ -16,8 +16,8 @@
*/
package org.orekit.rugged.core.duvenhage;
-import java.util.ArrayDeque;
-import java.util.Deque;
+import java.util.ArrayList;
+import java.util.List;
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.orekit.bodies.GeodeticPoint;
@@ -99,12 +99,12 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
// find where line-of-sight exit tile
final LimitPoint exit = findExit(tile, ellipsoid, position, los);
- final Deque<GeodeticPoint> lineOfSightQueue = new ArrayDeque<GeodeticPoint>();
- lineOfSightQueue.addFirst(exit.getPoint());
+ final List<GeodeticPoint> lineOfSightQueue = new ArrayList<GeodeticPoint>();
+ lineOfSightQueue.add(exit.getPoint());
while (!lineOfSightQueue.isEmpty()) {
- final GeodeticPoint next = lineOfSightQueue.removeFirst();
+ final GeodeticPoint next = lineOfSightQueue.remove(lineOfSightQueue.size() - 1);
final int nextLatIndex = tile.getLatitudeIndex(next.getLatitude());
final int nextLonIndex = tile.getLontitudeIndex(next.getLongitude());
if (currentLatIndex == nextLatIndex && currentLonIndex == nextLonIndex) {
@@ -123,11 +123,15 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
} else {
- // find the largest level in the min/max kd-tree at which entry and exit share a sub-tile
+ // 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) {
- // TODO: push intermediate points at sub-tiles boundaries on the queue
- throw RuggedException.createInternalError(null);
+ // 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,
+ nextLatIndex, nextLonIndex,
+ currentLatIndex, currentLonIndex));
} else {
if (next.getAltitude() >= tile.getMaxElevation(nextLatIndex, nextLonIndex, level)) {
// the line-of-sight segment is fully above Digital Elevation Model
@@ -139,30 +143,15 @@ public class DuvenhageAlgorithm implements IntersectionAlgorithm {
// search by using a finer-grained level in the min/max kd-tree
// push the point back into the queue
- lineOfSightQueue.addFirst(next);
-
- // increase level to have a more fine-grained min/max raster
- // (we go closer to individual pixels, or un-merge the merged min/max)
- // this un-merging implies finding the boundary between sub-tiles that
- // are merged at level l and not merged at level l+1. This boundary is
- // an iso-latitude if the merge is a row merging and is an iso-longitude
- // if the merge is a column merging. We therefore first identify the
- // row/column corresponding to the merging, then find the intersection of
- // the line-of-sight with the iso-surface, then insert this new intermediate
- // point in the queue
- final Vector3D mergePoint;
- if (tile.isColumnMerging(level + 1)) {
- final double longitudeMerge =
- tile.getLongitudeAtIndex(tile.getMergingColumn(currentLonIndex, level + 1));
- mergePoint = ellipsoid.pointAtLongitude(position, los, longitudeMerge);
- } else {
- final double latitudeSplit =
- tile.getLatitudeAtIndex(tile.getMergingRow(currentLatIndex, level + 1));
- mergePoint = ellipsoid.pointAtLatitude(position, los, latitudeSplit);
- }
-
- // push the intermediate point at sub-tile split into the queue
- lineOfSightQueue.addFirst(ellipsoid.transform(mergePoint, ellipsoid.getBodyFrame(), null));
+ 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,
+ nextLatIndex, nextLonIndex,
+ currentLatIndex, currentLonIndex));
// the current point remains the same
@@ -198,6 +187,56 @@ 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 boundary
+ * between sub-tiles that are merged at level l and not merged at level l+1.
+ * This boundary is an iso-latitude if the merge is a row merging and is an
+ * iso-longitude if the merge is a column merging.
+ * </p>
+ * @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 nextLatitude latitude index of next point (closer to Earth)
+ * @param nextLongitude longitude index of next point (closer to Earth)
+ * @param currentLatitude latitude index of current point (closer to satellite)
+ * @param currentLongitude longitude index of current point (closer to satellite)
+ * @return point corresponding to line-of-sight crossing the longitude/latitude
+ * limit between the un-merged sub-tiles at level-1
+ * @exception RuggedException if intersection point cannot be computed
+ * @exception OrekitException if intersection point 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 nextLatitude, final int nextLongitude,
+ final int currentLatitude, final int currentLongitude)
+ throws RuggedException, OrekitException {
+
+ final List<GeodeticPoint> crossings = new ArrayList<GeodeticPoint>();
+
+ if (tile.isColumnMerging(level + 1)) {
+ // sub-tiles at current level come from column merging at deeper level
+ for (final int longitudeIndex : tile.getCrossedBoundaryColumns(nextLongitude, currentLongitude, level)) {
+ final double crossingLongitude = tile.getLongitudeAtIndex(longitudeIndex);
+ final Vector3D crossingPoint = ellipsoid.pointAtLongitude(position, los, crossingLongitude);
+ crossings.add(ellipsoid.transform(crossingPoint, ellipsoid.getBodyFrame(), null));
+ }
+ } else {
+ // sub-tiles at current level come from row merging at deeper level
+ for (final int latitudeIndex : tile.getCrossedBoundaryRows(nextLatitude, currentLatitude, level)) {
+ final double crossingLatitude = tile.getLatitudeAtIndex(latitudeIndex);
+ final Vector3D crossingPoint = ellipsoid.pointAtLatitude(position, los, crossingLatitude);
+ crossings.add(ellipsoid.transform(crossingPoint, ellipsoid.getBodyFrame(), null));
+ }
+ }
+
+ return crossings;
+
+ }
+
/** Compute a line-of-sight exit point from a tile.
* @param tile tile to consider
* @param ellipsoid reference ellipsoid
diff --git a/rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTile.java b/rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTile.java
index 6317f76b..89e9f727 100644
--- a/rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTile.java
+++ b/rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTile.java
@@ -19,10 +19,57 @@ package org.orekit.rugged.core.duvenhage;
import org.apache.commons.math3.analysis.BivariateFunction;
import org.apache.commons.math3.analysis.function.Max;
import org.apache.commons.math3.analysis.function.Min;
+import org.apache.commons.math3.util.FastMath;
import org.orekit.rugged.core.raster.SimpleTile;
import org.orekit.rugged.core.raster.Tile;
/** Simple implementation of a {@link Tile} with a min/max kd tree.
+ * <p>
+ * A n level min/max kd-tree contains sub-tiles merging individual pixels
+ * together from coarse-grained (at level 0) to fine-grained (at level n-1).
+ * Level n-1, which is the deepest one, is computed from the raw pixels by
+ * merging adjacent pixels pairs columns (i.e. pixels at indices (i, 2j)
+ * and (i, 2j+1) are merged together by computing and storing the minimum
+ * and maxium in a sub-tile. Level n-1 therefore has the same number of rows
+ * but half the number of columns of the raw tile, and its sub-tiles are
+ * 1 pixel high and 2 pixels wide. Level n-2 is computed from level n-1 by
+ * merging sub-tiles rows. Level n-2 therefore has half the number of rows
+ * and half the number of columns of the raw tile, and its sub-tiles are
+ * 2 pixels high and 2 pixels wide. Level n-3 is again computed by merging
+ * columns, level n-4 merging rows and so on. As depth decreases, the number
+ * of sub-tiles decreases and their size increase. Level 0 is reached when
+ * there is only either one row or one column of large sub-tiles.
+ * </p>
+ * <p>
+ * During the merging process, if at some stage there is an odd number of
+ * rows or columns, then the last sub-tile at next level will not be computed
+ * by merging two rows/columns from the current level, but instead computed by
+ * simply copying the last single row/column. The process is therefore well
+ * defined for any raw tile initial dimensions. A direct consequence is that
+ * the dimension of the sub-tiles in the last row or column may be smaller than
+ * the dimension of regular sub-tiles.
+ * </p>
+ * <p>
+ * If we consider for example a tall 107 ⨉ 19 raw tile, the min/max kd-tree will
+ * have 9 levels:
+ * </p>
+ * <p>
+ * <table border="0">
+ * <tr BGCOLOR="#EEEEFF"><font size="+1">
+ * <td>Level</td> <td>Number of sub-tiles</td> <td>Regular sub-tiles dimension</td></font></tr>
+ * <tr> <td align="center">8</td> <td align="center">107 ⨉ 10</td> <td align="center"> 1 ⨉ 2</td>
+ * <tr> <td align="center">7</td> <td align="center"> 54 ⨉ 10</td> <td align="center"> 2 ⨉ 2</td>
+ * <tr> <td align="center">6</td> <td align="center"> 54 ⨉ 5</td> <td align="center"> 2 ⨉ 4</td>
+ * <tr> <td align="center">5</td> <td align="center"> 27 ⨉ 5</td> <td align="center"> 4 ⨉ 4</td>
+ * <tr> <td align="center">4</td> <td align="center"> 27 ⨉ 3</td> <td align="center"> 4 ⨉ 8</td>
+ * <tr> <td align="center">3</td> <td align="center"> 14 ⨉ 3</td> <td align="center"> 8 ⨉ 8</td>
+ * <tr> <td align="center">2</td> <td align="center"> 14 ⨉ 2</td> <td align="center"> 8 ⨉ 16</td>
+ * <tr> <td align="center">1</td> <td align="center"> 7 ⨉ 2</td> <td align="center">16 ⨉ 16</td>
+ * <tr> <td align="center">0</td> <td align="center"> 7 ⨉ 1</td> <td align="center">16 ⨉ 32</td>
+ * </table>
+ * </p>
+
+ * </p>
* @see MinMaxTreeTileFactory
* @author Luc Maisonobe
*/
@@ -121,13 +168,13 @@ public class MinMaxTreeTile extends SimpleTile {
}
- /** Get the largest level at which two pixels are merged in the same min/max sub-tile.
+ /** Get the deepest level at which two pixels are merged in the same min/max sub-tile.
* @param i1 row index of first pixel
* @param j1 column index of first pixel
* @param i2 row index of second pixel
* @param j2 column index of second pixel
- * @return largest level at which two pixels are merged in the same min/max sub-tile,
- * or negative if they are never merged in the same sub-tile
+ * @return deepest level at which two pixels are merged in the same min/max sub-tile,
+ * or -1 if they are never merged in the same sub-tile
* @see #getLevels()
* @see #getMinElevation(int, int, int)
* @see #getMaxElevation(int, int, int)
@@ -155,40 +202,76 @@ public class MinMaxTreeTile extends SimpleTile {
}
- /** Get the row at which two sub-tiles at level l are merged to give one sub-tile at level l-1.
+ /** Get the index of sub-tiles start rows crossed.
* <p>
- * This method is expected to be called for levels for which {@link #isColumnMerging(int)
- * isColumnMerging(level)} returns {@code false}. Calling it for other levels returns
- * unspecified results.
+ * When going from one row to another row at some tree level,
+ * we cross sub-tiles boundaries. This method returns the index
+ * of these boundaries.
* </p>
- * @param i row index of pixel in current sub-tile
- * @param level tree level to be merged into a lower level
- * @return index of row at which higher level sub-tiles were merged
- * (beware that this may be {@link #getLatitudeRows()} or more if the last row is not
- * really merged because only one pixel is available at this level)
+ * @param row1 starting row
+ * @param row2 ending row
+ * @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
*/
- public int getMergingRow(final int i, final int level) {
- final int k = start.length + 1 - level;
- final int rowShift = k / 2;
- return (i & (-1 << rowShift)) + (1 << (rowShift - 1));
+ public int[] getCrossedBoundaryRows(final int row1, final int row2, final int level) {
+
+ // number of rows in each sub-tile
+ final int rows = 1 << ((start.length - level) / 2);
+
+ if (row1 <= row2) {
+ return buildCrossings(rows * (row1 / rows + 1), row2, rows);
+ } else {
+ return buildCrossings(rows * ((row1 - 1) / rows), row2, -rows);
+ }
+
}
- /** Get the column at which two sub-tiles at level l are merged to give one sub-tile at level l-1.
+ /** Get the index of sub-tiles start columns crossed.
* <p>
- * This method is expected to be called for levels for which {@link #isColumnMerging(int)
- * isColumnMerging(level)} returns {@code true}. Calling it for other levels returns
- * unspecified results.
+ * When going from one column to another column at some tree level,
+ * we cross sub-tiles boundaries. This method returns the index
+ * of these boundaries.
* </p>
- * @param j column index of pixel in current sub-tile
- * @param level tree level to be merged into a lower level
- * @return index of column at which higher level sub-tiles were merged
- * (beware that this may be {@link #getLongitudeColumns()} or more if the last column is not
- * really merged because only one pixel is available at this level)
+ * @param column1 starting column
+ * @param column2 ending column
+ * @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 getMergingColumn(final int j, final int level) {
- final int k = start.length + 1 - level;
- final int colShift = (k + 1) / 2;
- return (j & (-1 << colShift)) + (1 << (colShift - 1));
+ public int[] getCrossedBoundaryColumns(final int column1, final int column2, final int level) {
+
+ // number of columns in each sub-tile
+ final int columns = 1 << ((start.length + 1 - level) / 2);;
+
+ if (column1 <= column2) {
+ return buildCrossings(columns * (column1 / columns + 1), column2, columns);
+ } else {
+ return buildCrossings(columns * ((column1 - 1) / columns), column2, -columns);
+ }
+
+ }
+
+ /** Build crossings arrays.
+ * @param begin begin crossing index
+ * @param end end crossing index (excluded, if equal to begin, the array is empty)
+ * @param step crossing step (may be negative)
+ * @return indices of rows or columns crossed at sub-tiles boundaries, in crossing order
+ */
+ private int[] buildCrossings(final int begin, final int end, final int step) {
+
+ // allocate array
+ final int[] crossings = new int[FastMath.max(0, (end - begin) / step)];
+
+ // fill it up
+ int crossing = begin;
+ for (int i = 0; i < crossings.length; ++i) {
+ crossings[i] = crossing;
+ crossing += step;
+ }
+
+ return crossings;
+
}
/** Check if the merging operation between level and level-1 is a column merging.
diff --git a/rugged-core/src/test/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java b/rugged-core/src/test/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java
index af10df10..d9afa570 100644
--- a/rugged-core/src/test/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java
+++ b/rugged-core/src/test/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java
@@ -62,10 +62,6 @@ public class MinMaxTreeTileTest {
maxTreeField.setAccessible(true);
Assert.assertEquals(2187, ((double[]) maxTreeField.get(tile)).length);
- for (int i = 0; i < tile.getLatitudeRows(); ++i) {
- Assert.assertEquals(32 * (i / 32) + 16, tile.getMergingRow(i, 0));
- }
-
}
@Test
@@ -189,43 +185,40 @@ public class MinMaxTreeTileTest {
}
@Test
- public void testMergingRow() throws RuggedException {
+ public void testSubTilesLimits() throws RuggedException {
for (int nbRows = 1; nbRows < 25; nbRows++) {
for (int nbColumns = 1; nbColumns < 25; nbColumns++) {
MinMaxTreeTile tile = createTile(nbRows, nbColumns);
- for (int i = 0; i < nbRows; i++) {
- for (int level = 0; level < tile.getLevels(); ++level) {
- if (!tile.isColumnMerging(level)) {
- int iMerge = tile.getMergingRow(i, level);
- if (iMerge < tile.getLatitudeRows()) {
- int[] neighbors1 = neighbors(iMerge - 1, 0, nbRows, nbColumns, tile.getLevels() - level);
- int[] neighbors2 = neighbors(iMerge, 0, nbRows, nbColumns, tile.getLevels() - level);
+ for (int level = 0; level < tile.getLevels(); ++level) {
+ for (int i1 = 0; i1 < nbRows; ++i1) {
+ for (int i2 = 0; i2 < nbRows; ++i2) {
+ int[] crossings = tile.getCrossedBoundaryRows(i1, i2, level);
+ if (FastMath.abs(i2 - i1) < 2) {
+ Assert.assertEquals(0, crossings.length);
+ }
+ for (int crossed : crossings) {
+ Assert.assertNotEquals(i1, crossed);
+ Assert.assertNotEquals(i2, crossed);
+ int[] neighbors1 = neighbors(crossed - 1, 0, nbRows, nbColumns, tile.getLevels() - level);
+ int[] neighbors2 = neighbors(crossed, 0, nbRows, nbColumns, tile.getLevels() - level);
Assert.assertEquals(neighbors1[1], neighbors2[0]);
}
}
- }
- }
- }
- }
- }
-
- @Test
- public void testMergingColumn() throws RuggedException {
- for (int nbRows = 1; nbRows < 25; nbRows++) {
- for (int nbColumns = 1; nbColumns < 25; nbColumns++) {
-
- MinMaxTreeTile tile = createTile(nbRows, nbColumns);
-
- for (int j = 0; j < nbColumns; j++) {
- for (int level = 0; level < tile.getLevels(); ++level) {
- if (tile.isColumnMerging(level)) {
- int jMerge = tile.getMergingColumn(j, level);
- if (jMerge < tile.getLongitudeColumns()) {
- int[] neighbors1 = neighbors(0, jMerge - 1, nbRows, nbColumns, tile.getLevels() - level);
- int[] neighbors2 = neighbors(0, jMerge, nbRows, nbColumns, tile.getLevels() - level);
- Assert.assertEquals(neighbors1[3], neighbors2[2]);
+ for (int j1 = 0; j1 < nbColumns; ++j1) {
+ for (int j2 = 0; j2 < nbColumns; ++j2) {
+ int[] crossings = tile.getCrossedBoundaryColumns(j1, j2, level);
+ if (FastMath.abs(j2 - j1) < 2) {
+ Assert.assertEquals(0, crossings.length);
+ }
+ for (int crossed : crossings) {
+ Assert.assertNotEquals(j1, crossed);
+ Assert.assertNotEquals(j2, crossed);
+ int[] neighbors1 = neighbors(0, crossed - 1, nbRows, nbColumns, tile.getLevels() - level);
+ int[] neighbors2 = neighbors(0, crossed, nbRows, nbColumns, tile.getLevels() - level);
+ Assert.assertEquals(neighbors1[3], neighbors2[2]);
+ }
}
}
}
--
GitLab