Added line-of-sight splitting at sub-tiles boundaries.

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

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.

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]);
+                                }
                             }
                         }
                     }