Started implementation of min/max KD-tree tile.

This is work in progress, the min/max computation seems wrong for now.

rugged-api/src/test/java/org/orekit/rugged/api/RuggedMessagesTest.java

@@ -29,7 +29,7 @@ public class RuggedMessagesTest {
 
     @Test
     public void testMessageNumber() {
-        Assert.assertEquals(3, RuggedMessages.values().length);
+        Assert.assertEquals(4, RuggedMessages.values().length);
     }
 
     @Test

rugged-core/src/main/java/org/orekit/rugged/core/dem/SimpleTile.java

@@ -87,6 +87,18 @@ public class SimpleTile implements Tile {
     /** {@inheritDoc} */
     @Override
     public void tileUpdateCompleted() throws RuggedException {
+        processUpdatedElevation(elevations);
+    }
+
+    /** Process elevation array at completion.
+     * </p>
+     * This method is called at tile update completion, it is
+     * expected to be overridden by subclasses. The default
+     * implementation does nothing.
+     * </p>
+     * @param elevations elevations array
+     */
+    protected void processUpdatedElevation(final double[] elevations) {
         // do nothing by default
     }
 

rugged-core/src/main/java/org/orekit/rugged/core/duvenhage/MinMaxTreeTile.java

@@ -16,7 +16,9 @@
  */
 package org.orekit.rugged.core.duvenhage;
 
-import org.orekit.rugged.api.RuggedException;
+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.orekit.rugged.core.dem.SimpleTile;
 import org.orekit.rugged.core.dem.Tile;
 
@@ -26,6 +28,21 @@ import org.orekit.rugged.core.dem.Tile;
  */
 public class MinMaxTreeTile extends SimpleTile {
 
+    /** Min kd-tree. */
+    private double[] minTree;
+
+    /** Max kd-tree. */
+    private double[] maxTree;
+
+    /** Start indices of tree levels. */
+    private int[] start;
+
+    /** Number of rows of tree levels. */
+    private int[] rows;
+
+    /** Number of columns of tree levels. */
+    private int[] columns;
+
     /** Simple constructor.
      * <p>
      * Creates an empty tile.
@@ -36,9 +53,201 @@ public class MinMaxTreeTile extends SimpleTile {
 
     /** {@inheritDoc} */
     @Override
-    public void tileUpdateCompleted() throws RuggedException {
-        // TODO: compute min/max tree
-        throw RuggedException.createInternalError(null);
+    protected void processUpdatedElevation(final double[] elevations) {
+
+        // set up the levels
+        final int size = setLevels(0, getLatitudeRows(), getLongitudeColumns());
+        minTree = new double[size];
+        maxTree = new double[size];
+
+        // compute min/max trees
+        applyRecursively(minTree, 0, getLatitudeRows(), getLongitudeColumns(),
+                         new Min(), elevations, 0);
+        applyRecursively(maxTree, 0, getLatitudeRows(), getLongitudeColumns(),
+                         new Max(), elevations, 0);
+
+    }
+
+    /** Get the number of kd-tree levels (not counting raw elevations).
+     * @return number of kd-tree levels
+     * @see #getMinElevation(int, int, int)
+     * @see #getMaxElevation(int, int, int)
+     */
+    public int getLevels() {
+        return start.length;
+    }
+
+    /** Get the minimum elevation at some level tree.
+     * @param i row index of the pixel
+     * @param j column index of the pixel
+     * @param level tree level
+     * @return minimum elevation
+     * @see #getLevels()
+     * @see #getMaxElevation(int, int, int)
+     */
+    public double getMinElevation(final int i, final int j, final int level) {
+
+        // compute row index in level merged array
+        final int levelI = i >> ((level + 1) / 2);
+        final int levelJ = j >> ((level + 2) / 2);
+
+        return minTree[start[level] + levelI * columns[level] + levelJ];
+
+    }
+
+    /** Get the maximum elevation at some level tree.
+     * @param i row index of the pixel
+     * @param j column index of the pixel
+     * @param level tree level
+     * @return maximum elevation
+     * @see #getLevels()
+     * @see #getMinElevation(int, int, int)
+     */
+    public double getMaxElevation(final int i, final int j, final int level) {
+
+        // compute row index in level merged array
+        final int levelI = i >> ((level + 1) / 2);
+        final int levelJ = j >> (level / 2);
+
+        return maxTree[start[level] + levelI * columns[level] + levelJ];
+
+    }
+
+    /** Recursive setting of tree levels.
+     * <p>
+     * The following algorithms works for any array shape, even with
+     * rows or columns which are not powers of 2 or with one
+     * dimension much larger than the other. As an example, starting
+     * from a 107 ⨉ 19 array, we get the following 9 levels, for a
+     * total of 2187 elements in both trees:
+     * </p>
+     * <p>
+     * <table border="0">
+     * <tr BGCOLOR="#EEEEFF"><font size="+1">
+     *     <td>Level</td>   <td>Dimension</td>  <td>Start index</td>  <td>End index</td></font></tr>
+     * <tr>   <td>8</td>     <td>  7 ⨉  1</td>       <td>   0</td>        <td>  6</td> </tr>
+     * <tr>   <td>7</td>     <td>  7 ⨉  2</td>       <td>   7</td>        <td> 20</td> </tr>
+     * <tr>   <td>6</td>     <td> 14 ⨉  2</td>       <td>  21</td>        <td> 48</td> </tr>
+     * <tr>   <td>5</td>     <td> 14 ⨉  3</td>       <td>  49</td>        <td> 90</td> </tr>
+     * <tr>   <td>4</td>     <td> 27 ⨉  3</td>       <td>  91</td>        <td>171</td> </tr>
+     * <tr>   <td>3</td>     <td> 27 ⨉  5</td>       <td> 172</td>        <td>306</td> </tr>
+     * <tr>   <td>2</td>     <td> 54 ⨉  5</td>       <td> 307</td>        <td>576</td> </tr>
+     * <tr>   <td>1</td>     <td> 54 ⨉ 10</td>      <td> 577</td>        <td>1116</td> </tr>
+     * <tr>   <td>0</td>     <td>107 ⨉ 10</td>      <td>1117</td>        <td>2186</td> </tr>
+     * </table>
+     * </p>
+     * @param level current level (counting from leafs to root)
+     * @param levelRows number of rows at current level
+     * @param levelColumns number of columns at current level
+     * @return size cumulative size from current level to root
+     */
+    private int setLevels(final int level, final int levelRows, final int levelColumns) {
+
+        if (levelRows == 1 || levelColumns == 1) {
+            // we have found root, stop recursion
+            start              = new int[level];
+            rows               = new int[level];
+            columns            = new int[level];
+            start[level - 1]   = 0;
+            rows[level - 1]    = levelRows;
+            columns[level - 1] = levelColumns;
+            return levelRows * levelColumns;
+        }
+
+        final int size;
+        if ((level & 0x1) == 0) {
+            // columns merging
+            size = setLevels(level + 1, levelRows, (levelColumns + 1) / 2);
+        } else {
+            // rows merging
+            size = setLevels(level + 1, (levelRows + 1) / 2, levelColumns);
+        }
+
+        if (level > 0) {
+            // store current level characteristics
+            start[level - 1]   = size;
+            rows[level - 1]    = levelRows;
+            columns[level - 1] = levelColumns;
+            return size + levelRows * levelColumns;
+        } else {
+            // we don't count the elements at leaf as they are not stored
+            // in the min/max trees
+            return size;
+        }
+
+    }
+
+    /** Recursive application of a function.
+     * @param tree to fill-up with the recursive applications
+     * @param level current level
+     * @param levelRows number of rows at current level
+     * @param levelColumns number of columns at current level
+     * @param f function to apply
+     * @param base base array from which function arguments are drawn
+     * @param first index of the first element to consider in base array
+     */
+    private void applyRecursively(final double[] tree,
+                                  final int level, final int levelRows, final int levelColumns,
+                                  final BivariateFunction f,
+                                  final double[] base, final int first) {
+        if ((level & 0x1) != (start.length & 0x1)) {
+
+            // merge columns pairs
+            int           iTree       = start[level];
+            int           iBase       = first;
+            final int     nextColumns = (levelColumns + 1) / 2;
+            final boolean odd         = (levelColumns & 0x1) != 0;
+            for (int i = 0; i < levelRows; ++i) {
+
+                // regular pairs
+                int jEnd = odd ? nextColumns - 1 : nextColumns;
+                for (int j = 0; j < jEnd; ++j) {
+                    tree[iTree++] = f.value(base[iBase], base[iBase + 1]);
+                    iBase += 2;
+                }
+
+                if (odd) {
+                    // last column
+                    tree[iTree++] = base[iBase++];
+                }
+
+
+            }
+
+            if (level < start.length - 1) {
+                applyRecursively(tree, level + 1, levelRows, nextColumns, f, tree, start[level]);
+            }
+
+        } else {
+
+            // merge rows pairs
+            int           iTree    = start[level];
+            int           iBase    = first;
+            final int     nextRows = (levelRows + 1) / 2;
+            final boolean odd      = (levelRows & 0x1) != 0;
+
+            // regular pairs
+            int iEnd = odd ? nextRows - 1 : nextRows;
+            for (int i = 0; i < iEnd; ++i) {
+
+                for (int j = 0; j < levelColumns; ++j) {
+                    tree[iTree++] = f.value(base[iBase], base[iBase + levelColumns]);
+                    iBase++;
+                }
+                iBase += levelColumns;
+
+            }
+
+            if (odd) {
+                // last row
+                System.arraycopy(base, iBase, tree, iTree, levelColumns);
+            }
+
+            if (level < start.length - 1) {
+                applyRecursively(tree, level + 1, nextRows, levelColumns, f, tree, start[level]);
+            }
+
+        }
     }
 
 }

rugged-core/src/test/java/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java

+/* Copyright 2013-2014 CS Systèmes d'Information
+ * Licensed to CS Systèmes d'Information (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package orekit.rugged.core.duvenhage;
+
+import java.lang.reflect.Field;
+
+import org.apache.commons.math3.random.RandomGenerator;
+import org.apache.commons.math3.random.Well19937a;
+import org.apache.commons.math3.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+import org.orekit.rugged.api.RuggedException;
+import org.orekit.rugged.core.duvenhage.MinMaxTreeTile;
+import org.orekit.rugged.core.duvenhage.MinMaxTreeTileFactory;
+
+public class MinMaxTreeTileTest {
+
+    @Test
+    public void testSizeTall()
+        throws RuggedException, SecurityException, NoSuchFieldException,
+               IllegalArgumentException, IllegalAccessException {
+        MinMaxTreeTile tile = new MinMaxTreeTileFactory().createTile();
+        tile.setGeometry(1.0, 2.0, 0.1, 0.2, 107, 19);
+        tile.tileUpdateCompleted();
+        Assert.assertEquals(9, tile.getLevels());
+
+        Field startField = MinMaxTreeTile.class.getDeclaredField("start");
+        startField.setAccessible(true);
+        int[] start = (int[]) startField.get(tile);
+        Assert.assertEquals(   0, start[ 8]);
+        Assert.assertEquals(   7, start[ 7]);
+        Assert.assertEquals(  21, start[ 6]);
+        Assert.assertEquals(  49, start[ 5]);
+        Assert.assertEquals(  91, start[ 4]);
+        Assert.assertEquals( 172, start[ 3]);
+        Assert.assertEquals( 307, start[ 2]);
+        Assert.assertEquals( 577, start[ 1]);
+        Assert.assertEquals(1117, start[ 0]);
+
+        Field rowsField = MinMaxTreeTile.class.getDeclaredField("rows");
+        rowsField.setAccessible(true);
+        int[] rows = (int[]) rowsField.get(tile);
+        Assert.assertEquals(   7, rows[ 8]);
+        Assert.assertEquals(   7, rows[ 7]);
+        Assert.assertEquals(  14, rows[ 6]);
+        Assert.assertEquals(  14, rows[ 5]);
+        Assert.assertEquals(  27, rows[ 4]);
+        Assert.assertEquals(  27, rows[ 3]);
+        Assert.assertEquals(  54, rows[ 2]);
+        Assert.assertEquals(  54, rows[ 1]);
+        Assert.assertEquals( 107, rows[ 0]);
+
+        Field columnsField = MinMaxTreeTile.class.getDeclaredField("columns");
+        columnsField.setAccessible(true);
+        int[] columns = (int[]) columnsField.get(tile);
+        Assert.assertEquals(  1, columns[ 8]);
+        Assert.assertEquals(  2, columns[ 7]);
+        Assert.assertEquals(  2, columns[ 6]);
+        Assert.assertEquals(  3, columns[ 5]);
+        Assert.assertEquals(  3, columns[ 4]);
+        Assert.assertEquals(  5, columns[ 3]);
+        Assert.assertEquals(  5, columns[ 2]);
+        Assert.assertEquals( 10, columns[ 1]);
+        Assert.assertEquals( 10, columns[ 0]);
+        Field minTreeField = MinMaxTreeTile.class.getDeclaredField("minTree");
+        minTreeField.setAccessible(true);
+        Assert.assertEquals(2187, ((double[]) minTreeField.get(tile)).length);
+        Field maxTreeField = MinMaxTreeTile.class.getDeclaredField("maxTree");
+        maxTreeField.setAccessible(true);
+        Assert.assertEquals(2187, ((double[]) maxTreeField.get(tile)).length);
+    }
+
+    @Test
+    public void testSizeFat()
+        throws RuggedException, SecurityException, NoSuchFieldException,
+               IllegalArgumentException, IllegalAccessException {
+        MinMaxTreeTile tile = new MinMaxTreeTileFactory().createTile();
+        tile.setGeometry(1.0, 2.0, 0.1, 0.2, 4, 7);
+        tile.tileUpdateCompleted();
+        Assert.assertEquals(4, tile.getLevels());
+
+        Field startField = MinMaxTreeTile.class.getDeclaredField("start");
+        startField.setAccessible(true);
+        int[] start = (int[]) startField.get(tile);
+        Assert.assertEquals( 0, start[ 3]);
+        Assert.assertEquals( 2, start[ 2]);
+        Assert.assertEquals( 6, start[ 1]);
+        Assert.assertEquals(14, start[ 0]);
+
+        Field rowsField = MinMaxTreeTile.class.getDeclaredField("rows");
+        rowsField.setAccessible(true);
+        int[] rows = (int[]) rowsField.get(tile);
+        Assert.assertEquals( 1, rows[ 3]);
+        Assert.assertEquals( 2, rows[ 2]);
+        Assert.assertEquals( 2, rows[ 1]);
+        Assert.assertEquals( 4, rows[ 0]);
+
+        Field columnsField = MinMaxTreeTile.class.getDeclaredField("columns");
+        columnsField.setAccessible(true);
+        int[] columns = (int[]) columnsField.get(tile);
+        Assert.assertEquals( 2, columns[ 3]);
+        Assert.assertEquals( 2, columns[ 2]);
+        Assert.assertEquals( 4, columns[ 1]);
+        Assert.assertEquals( 4, columns[ 0]);
+
+        Field minTreeField = MinMaxTreeTile.class.getDeclaredField("minTree");
+        minTreeField.setAccessible(true);
+        Assert.assertEquals(30, ((double[]) minTreeField.get(tile)).length);
+        Field maxTreeField = MinMaxTreeTile.class.getDeclaredField("maxTree");
+        maxTreeField.setAccessible(true);
+        Assert.assertEquals(30, ((double[]) maxTreeField.get(tile)).length);
+    }
+
+    @Test
+    public void testMinMax() throws RuggedException {
+        RandomGenerator random = new Well19937a(0xfbbc1d1739b23555l);
+        checkMinMax(4, 7, 100, random);
+    }
+
+    private void checkMinMax(int nbRows, int nbColumns, int nbChecks, RandomGenerator random)
+        throws RuggedException {
+
+        MinMaxTreeTile tile = new MinMaxTreeTileFactory().createTile();
+        tile.setGeometry(1.0, 2.0, 0.1, 0.2, nbRows, nbColumns);
+        for (int i = 0; i < nbRows; ++i) {
+            for (int j = 0; j < nbColumns; ++j) {
+                tile.setElevation(i, j, 1000 * random.nextDouble());
+            }
+        }
+        tile.tileUpdateCompleted();
+
+        for (int k = 0; k < nbChecks; ++k) {
+            int row    = random.nextInt(nbRows);
+            int column = random.nextInt(nbColumns);
+            int level  = random.nextInt(tile.getLevels());
+
+            // reference min and max
+            int[] neighbors = neighbors(row, column, nbRows, nbColumns, level);
+            System.out.println(row + " " + column + " (" + level + "): [" +
+                               neighbors[0] + ", " + neighbors[1] + "] [" +
+                               neighbors[2] + ", " + neighbors[3] + "]");
+            double min = Double.POSITIVE_INFINITY;
+            double max = Double.POSITIVE_INFINITY;
+            for (int i = neighbors[0]; i < neighbors[1]; ++i) {
+                for (int j = neighbors[2]; j < neighbors[3]; ++j) {
+                    double pixelValue = tile.getElevationAtIndices(i, j);
+                    min = FastMath.min(min, pixelValue);
+                    max = FastMath.max(max, pixelValue);
+                }
+            }
+
+            Assert.assertEquals(min, tile.getMinElevation(row, column, level), 1.0e-10 * min);
+            Assert.assertEquals(max, tile.getMaxElevation(row, column, level), 1.0e-10 * max);
+
+        }
+
+    }
+
+    private int[] neighbors(int row, int column, int nbRows, int nbColumns, int level) {
+
+        // poor man identification of neighbors cells merged together with specified cell
+        int rMin = 0;
+        int rMax = row;
+        int cMin = 0;
+        int cMax = column;
+
+        boolean mergeColumns = true;
+        for (int i = 0; i <= level; ++i) {
+            if (mergeColumns) {
+                int split = (nbColumns + 1) / 2;
+                if (column < split) {
+                    cMax = split;
+                } else {
+                    cMin = split;
+                }
+                nbColumns = cMax - cMin;
+            } else {
+                int split = (nbRows + 1) / 2;
+                if (row < split) {
+                    rMax = split;
+                } else {
+                    rMin = split;
+                }
+                nbRows = rMax - rMin;
+            }
+            mergeColumns = !mergeColumns;
+        }
+
+        return new int[] { rMin, rMax, cMin, cMax };
+
+    }
+
+}