From 5940ea11acdb3f4fe79223e803f0ff395ea0b00e Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Mon, 10 Mar 2014 18:38:08 +0100
Subject: [PATCH] Started implementation of min/max KD-tree tile.
This is work in progress, the min/max computation seems wrong for now.
---
.../orekit/rugged/api/RuggedMessagesTest.java | 2 +-
.../orekit/rugged/core/dem/SimpleTile.java | 12 +
.../rugged/core/duvenhage/MinMaxTreeTile.java | 217 +++++++++++++++++-
.../core/duvenhage/MinMaxTreeTileTest.java | 207 +++++++++++++++++
4 files changed, 433 insertions(+), 5 deletions(-)
create mode 100644 rugged-core/src/test/java/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java
diff --git a/rugged-api/src/test/java/org/orekit/rugged/api/RuggedMessagesTest.java b/rugged-api/src/test/java/org/orekit/rugged/api/RuggedMessagesTest.java
index 543ac50c..b220c71b 100644
--- a/rugged-api/src/test/java/org/orekit/rugged/api/RuggedMessagesTest.java
+++ b/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
diff --git a/rugged-core/src/main/java/org/orekit/rugged/core/dem/SimpleTile.java b/rugged-core/src/main/java/org/orekit/rugged/core/dem/SimpleTile.java
index a6b6ed61..7a84e885 100644
--- a/rugged-core/src/main/java/org/orekit/rugged/core/dem/SimpleTile.java
+++ b/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
}
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 480a85a8..3dcf0b26 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
@@ -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]);
+ }
+
+ }
}
}
diff --git a/rugged-core/src/test/java/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java b/rugged-core/src/test/java/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java
new file mode 100644
index 00000000..68295ff6
--- /dev/null
+++ b/rugged-core/src/test/java/orekit/rugged/core/duvenhage/MinMaxTreeTileTest.java
@@ -0,0 +1,207 @@
+/* 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 };
+
+ }
+
+}
--
GitLab