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 543ac50cd16897f4c04ff83b5286530e072c0526..b220c71b2b5aba28d56fd71f7e78f1800210d8c7 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 a6b6ed61c33da44cfdda26f07193f68569627fb3..7a84e885e9468f9b417a16014cb63044acc5fa06 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 480a85a8a6087cf408b28b2e14f844a3a85cb870..3dcf0b26da4c44fd830602897f1509e9bbef7181 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 0000000000000000000000000000000000000000..68295ff633cd16a6f3e0c0aba4ea892f372fb9bf --- /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 }; + + } + +}