MinMaxTreeTile.java

  1. /* Copyright 2013-2022 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.rugged.intersection.duvenhage;

  19. import org.hipparchus.util.FastMath;
  20. import org.orekit.rugged.errors.DumpManager;
  21. import org.orekit.rugged.raster.SimpleTile;
  22. import org.orekit.rugged.utils.MaxSelector;
  23. import org.orekit.rugged.utils.MinSelector;
  24. import org.orekit.rugged.utils.Selector;

  26. /** Implementation of a {@link org.orekit.rugged.raster.Tile} with a min/max kd tree.
  27.  * <p>
  28.  * A n level min/max kd-tree contains sub-tiles merging individual cells
  29.  * together from coarse-grained (at level 0) to fine-grained (at level n-1).
  30.  * Level n-1, which is the deepest one, is computed from the raw cells by
  31.  * merging adjacent cells pairs columns (i.e. cells at indices (i, 2j)
  32.  * and (i, 2j+1) are merged together by computing and storing the minimum
  33.  * and maximum in a sub-tile. Level n-1 therefore has the same number of rows
  34.  * but half the number of columns of the raw tile, and its sub-tiles are
  35.  * 1 cell high and 2 cells wide. Level n-2 is computed from level n-1 by
  36.  * merging sub-tiles rows. Level n-2 therefore has half the number of rows
  37.  * and half the number of columns of the raw tile, and its sub-tiles are
  38.  * 2 cells high and 2 cells wide. Level n-3 is again computed by merging
  39.  * columns, level n-4 merging rows and so on. As depth decreases, the number
  40.  * of sub-tiles decreases and their size increase. Level 0 is reached when
  41.  * there is only either one row or one column of large sub-tiles.
  42.  * </p>
  43.  * <p>
  44.  * During the merging process, if at some stage there is an odd number of
  45.  * rows or columns, then the last sub-tile at next level will not be computed
  46.  * by merging two rows/columns from the current level, but instead computed by
  47.  * simply copying the last single row/column. The process is therefore well
  48.  * defined for any raw tile initial dimensions. A direct consequence is that
  49.  * the dimension of the sub-tiles in the last row or column may be smaller than
  50.  * the dimension of regular sub-tiles.
  51.  * </p>
  52.  * <p>
  53.  * If we consider for example a tall 107 ⨉ 19 raw tile, the min/max kd-tree will
  54.  * have 9 levels:
  55.  * </p>
  56.  *
  57.  * <table summary="" border="0">
  58.  * <tr style="background-color:#EEEEFF;">
  59.  *             <td>Level</td>         <td>Number of sub-tiles</td>    <td>Regular sub-tiles dimension</td></tr>
  60.  * <tr>  <td align="center">8</td>  <td align="center">107 ⨉ 10</td>       <td align="center"> 1 ⨉   2</td>
  61.  * <tr>  <td align="center">7</td>  <td align="center"> 54 ⨉ 10</td>       <td align="center"> 2 ⨉   2</td>
  62.  * <tr>  <td align="center">6</td>  <td align="center"> 54 ⨉  5</td>        <td align="center"> 2 ⨉  4</td>
  63.  * <tr>  <td align="center">5</td>  <td align="center"> 27 ⨉  5</td>        <td align="center"> 4 ⨉  4</td>
  64.  * <tr>  <td align="center">4</td>  <td align="center"> 27 ⨉  3</td>        <td align="center"> 4 ⨉  8</td>
  65.  * <tr>  <td align="center">3</td>  <td align="center"> 14 ⨉  3</td>        <td align="center"> 8 ⨉  8</td>
  66.  * <tr>  <td align="center">2</td>  <td align="center"> 14 ⨉  2</td>        <td align="center"> 8 ⨉ 16</td>
  67.  * <tr>  <td align="center">1</td>  <td align="center">  7 ⨉  2</td>        <td align="center">16 ⨉ 16</td>
  68.  * <tr>  <td align="center">0</td>  <td align="center">  7 ⨉  1</td>        <td align="center">16 ⨉ 32</td>
  69.  * </table>
  70.  * <p>
  71.  * @see MinMaxTreeTileFactory
  72.  * @author Luc Maisonobe
  73.  */
  74. public class MinMaxTreeTile extends SimpleTile {

  76.     /** Raw elevations. */
  77.     private double[] raw;

  79.     /** Min kd-tree. */
  80.     private double[] minTree;

  82.     /** Max kd-tree. */
  83.     private double[] maxTree;

  85.     /** Start indices of tree levels. */
  86.     private int[] start;

  88.     /** Simple constructor.
  89.      * <p>
  90.      * Creates an empty tile.
  91.      * </p>
  92.      */
  93.     MinMaxTreeTile() {
  94.     }

  96.     /** {@inheritDoc} */
  97.     @Override
  98.     protected void processUpdatedElevation(final double[] elevations) {

  100.         raw = elevations;

  102.         final int nbRows = getLatitudeRows();
  103.         final int nbCols = getLongitudeColumns();

  105.         // set up the levels
  106.         final int size = setLevels(0, nbRows, nbCols);
  107.         minTree = new double[size];
  108.         maxTree = new double[size];

  110.         // compute min/max trees
  111.         if (start.length > 0) {

  113.             final double[] preprocessed = new double[raw.length];

  115.             preprocess(preprocessed, raw, nbRows, nbCols, MinSelector.getInstance());
  116.             applyRecursively(minTree, start.length - 1, nbRows, nbCols, MinSelector.getInstance(), preprocessed, 0);

  118.             preprocess(preprocessed, raw, nbRows, nbCols, MaxSelector.getInstance());
  119.             applyRecursively(maxTree, start.length - 1, nbRows, nbCols, MaxSelector.getInstance(), preprocessed, 0);

  121.         }

  123.     }

  125.     /** Get the number of kd-tree levels (not counting raw elevations).
  126.      * @return number of kd-tree levels
  127.      * @see #getMinElevation(int, int, int)
  128.      * @see #getMaxElevation(int, int, int)
  129.      * @see #getMergeLevel(int, int, int, int)
  130.      */
  131.     public int getLevels() {
  132.         return start.length;
  133.     }

  135.     /** Get the minimum elevation at some level tree.
  136.      * <p>
  137.      * Note that the min elevation is <em>not</em> computed
  138.      * only at cell center, but considering that it is interpolated
  139.      * considering also Eastwards and Northwards neighbors, and extends
  140.      * up to the center of these neighbors. As an example, lets consider
  141.      * four neighboring cells in some Digital Elevation Model:
  142.      * <table summary="" border="0" cellpadding="5" style="background-color:#f5f5dc;">
  143.      * <tr><th style="background-color:#c9d5c9;">j+1</th><td>11</td><td>10</td></tr>
  144.      * <tr><th style="background-color:#c9d5c9;">j</th><td>12</td><td>11</td></tr>
  145.      * <tr  style="background-color:#c9d5c9;"><th>j/i</th><th>i</th><th>i+1</th></tr>
  146.      * </table>
  147.      * When we interpolate elevation at a point located slightly South-West
  148.      * to the center of the (i+1, j+1) cell, we use all four cells in the
  149.      * interpolation, and we will get a result very close to 10 if we start
  150.      * close to (i+1, j+1) cell center. As the min value for this interpolation
  151.      * is stored at (i, j) indices, this implies that {@code getMinElevation(i,
  152.      * j, l)} must return 10 if l is chosen such that the sub-tile at
  153.      * tree level l includes cell (i,j) but not cell (i+1, j+1). In other words,
  154.      * interpolation implies sub-tile boundaries are overshoot by one column to
  155.      * the East and one row to the North when computing min.
  156.      *
  157.      * @param i row index of the cell
  158.      * @param j column index of the cell
  159.      * @param level tree level
  160.      * @return minimum value that can be reached when interpolating elevation
  161.      * in the sub-tile
  162.      * @see #getLevels()
  163.      * @see #getMaxElevation(int, int, int)
  164.      * @see #getMergeLevel(int, int, int, int)
  165.      */
  166.     public double getMinElevation(final int i, final int j, final int level) {

  168.         // compute indices in level merged array
  169.         final int k        = start.length - level;
  170.         final int rowShift = k / 2;
  171.         final int colShift = (k + 1) / 2;
  172.         final int levelI   = i >> rowShift;
  173.         final int levelJ   = j >> colShift;
  174.         final int levelC   = 1 + ((getLongitudeColumns() - 1) >> colShift);

  176.         if (DumpManager.isActive()) {
  177.             final int[] min = locateMin(i, j, level);
  178.             final int index = min[0] * getLongitudeColumns() + min[1];
  179.             DumpManager.dumpTileCell(this, min[0],     min[1],     raw[index]);
  180.             if (index + getLongitudeColumns() < raw.length) {
  181.                 DumpManager.dumpTileCell(this, min[0] + 1, min[1],     raw[index + getLongitudeColumns()]);
  182.             }
  183.             if (index + 1 < raw.length) {
  184.                 DumpManager.dumpTileCell(this, min[0],     min[1] + 1, raw[index + 1]);
  185.             }
  186.             if (index + getLongitudeColumns() + 1 < raw.length) {
  187.                 DumpManager.dumpTileCell(this, min[0] + 1, min[1] + 1, raw[index + getLongitudeColumns() + 1]);
  188.             }
  189.         }

  191.         return minTree[start[level] + levelI * levelC + levelJ];

  193.     }

  195.     /** Get the maximum elevation at some level tree.
  196.      * <p>
  197.      * Note that the max elevation is <em>not</em> computed
  198.      * only at cell center, but considering that it is interpolated
  199.      * considering also Eastwards and Northwards neighbors, and extends
  200.      * up to the center of these neighbors. As an example, lets consider
  201.      * four neighboring cells in some Digital Elevation Model:
  202.      * <table summary="" border="0" cellpadding="5" style="background-color:#f5f5dc;">
  203.      * <tr><th style="background-color:#c9d5c9;">j+1</th><td>11</td><td>12</td></tr>
  204.      * <tr><th style="background-color:#c9d5c9;">j</th><td>10</td><td>11</td></tr>
  205.      * <tr  style="background-color:#c9d5c9;"><th>j/i</th><th>i</th><th>i+1</th></tr>
  206.      * </table>
  207.      * When we interpolate elevation at a point located slightly South-West
  208.      * to the center of the (i+1, j+1) cell, we use all four cells in the
  209.      * interpolation, and we will get a result very close to 12 if we start
  210.      * close to (i+1, j+1) cell center. As the max value for this interpolation
  211.      * is stored at (i, j) indices, this implies that {@code getMaxElevation(i,
  212.      * j, l)} must return 12 if l is chosen such that the sub-tile at
  213.      * tree level l includes cell (i,j) but not cell (i+1, j+1). In other words,
  214.      * interpolation implies sub-tile boundaries are overshoot by one column to
  215.      * the East and one row to the North when computing max.
  216.      *
  217.      * @param i row index of the cell
  218.      * @param j column index of the cell
  219.      * @param level tree level
  220.      * @return maximum value that can be reached when interpolating elevation
  221.      * in the sub-tile
  222.      * @see #getLevels()
  223.      * @see #getMinElevation(int, int, int)
  224.      * @see #getMergeLevel(int, int, int, int)
  225.      */
  226.     public double getMaxElevation(final int i, final int j, final int level) {

  228.         // compute indices in level merged array
  229.         final int k        = start.length - level;
  230.         final int rowShift = k / 2;
  231.         final int colShift = (k + 1) / 2;
  232.         final int levelI   = i >> rowShift;
  233.         final int levelJ   = j >> colShift;
  234.         final int levelC   = 1 + ((getLongitudeColumns() - 1) >> colShift);

  236.         if (DumpManager.isActive()) {
  237.             final int[] max = locateMax(i, j, level);
  238.             final int index = max[0] * getLongitudeColumns() + max[1];
  239.             DumpManager.dumpTileCell(this, max[0],     max[1],     raw[index]);
  240.             if (index + getLongitudeColumns() < raw.length) {
  241.                 DumpManager.dumpTileCell(this, max[0] + 1, max[1],     raw[index + getLongitudeColumns()]);
  242.             }
  243.             if (index + 1 < raw.length) {
  244.                 DumpManager.dumpTileCell(this, max[0],     max[1] + 1, raw[index + 1]);
  245.             }
  246.             if (index + getLongitudeColumns() + 1 < raw.length) {
  247.                 DumpManager.dumpTileCell(this, max[0] + 1, max[1] + 1, raw[index + getLongitudeColumns() + 1]);
  248.             }
  249.         }

  251.         return maxTree[start[level] + levelI * levelC + levelJ];

  253.     }

  255.     /** Locate the cell at which min elevation is reached for a specified level.
  256.      * <p>
  257.      * Min is computed with respect to the continuous interpolated elevation, which
  258.      * takes four neighboring cells into account. This implies that the cell at which
  259.      * min value is reached for some level is either within the sub-tile for this level,
  260.      * or in some case it may be one column outside to the East or one row outside to
  261.      * the North. See {@link #getMinElevation()} for a more complete explanation.
  262.      * </p>
  263.      * @param i row index of the cell
  264.      * @param j column index of the cell
  265.      * @param level tree level of the sub-tile considered
  266.      * @return row/column indices of the cell at which min elevation is reached
  267.      */
  268.     public int[] locateMin(final int i, final int j, final int level) {
  269.         return locateMinMax(i, j, level, MinSelector.getInstance(), minTree);
  270.     }

  272.     /** Locate the cell at which max elevation is reached for a specified level.
  273.      * <p>
  274.      * Max is computed with respect to the continuous interpolated elevation, which
  275.      * takes four neighboring cells into account. This implies that the cell at which
  276.      * max value is reached for some level is either within the sub-tile for this level,
  277.      * or in some case it may be one column outside to the East or one row outside to
  278.      * the North. See {@link #getMaxElevation()} for a more complete explanation.
  279.      * </p>
  280.      * @param i row index of the cell
  281.      * @param j column index of the cell
  282.      * @param level tree level of the sub-tile considered
  283.      * @return row/column indices of the cell at which min elevation is reached
  284.      */
  285.     public int[] locateMax(final int i, final int j, final int level) {
  286.         return locateMinMax(i, j, level, MaxSelector.getInstance(), maxTree);
  287.     }

  289.     /** Locate the cell at which min/max elevation is reached for a specified level.
  290.      * @param i row index of the cell
  291.      * @param j column index of the cell
  292.      * @param level tree level of the sub-tile considered
  293.      * @param selector min/max selector to use
  294.      * @param tree min/max tree to use
  295.      * @return row/column indices of the cell at which min/max elevation is reached
  296.      */
  297.     private int[] locateMinMax(final int i, final int j, final int level,
  298.                                final Selector selector, final double[] tree) {

  300.         final int k  = start.length - level;
  301.         int rowShift = k / 2;
  302.         int colShift = (k + 1) / 2;
  303.         int levelI   = i >> rowShift;
  304.         int levelJ   = j >> colShift;
  305.         int levelR   = 1 + ((getLatitudeRows()     - 1) >> rowShift);
  306.         int levelC   = 1 + ((getLongitudeColumns() - 1) >> colShift);

  308.         // track the cell ancestors from merged tree at specified level up to tree at level 1
  309.         for (int l = level + 1; l < start.length; ++l) {

  311.             if (isColumnMerging(l)) {

  313.                 --colShift;
  314.                 levelC = 1 + ((getLongitudeColumns() - 1) >> colShift);
  315.                 levelJ = levelJ << 1;

  317.                 if (levelJ + 1 < levelC) {
  318.                     // the cell results from a regular merging of two columns
  319.                     if (selector.selectFirst(tree[start[l] + levelI * levelC + levelJ + 1],
  320.                                              tree[start[l] + levelI * levelC + levelJ])) {
  321.                         levelJ++;
  322.                     }
  323.                 }

  325.             } else {

  327.                 --rowShift;
  328.                 levelR = 1 + ((getLatitudeRows() - 1) >> rowShift);
  329.                 levelI = levelI << 1;

  331.                 if (levelI + 1 < levelR) {
  332.                     // the cell results from a regular merging of two rows
  333.                     if (selector.selectFirst(tree[start[l] + (levelI + 1) * levelC + levelJ],
  334.                                              tree[start[l] + levelI       * levelC + levelJ])) {
  335.                         levelI++;
  336.                     }
  337.                 }

  339.             }

  341.         }

  343.         // we are now at first merge level, which always results from a column merge
  344.         // or pre-processed data, which themselves result from merging four cells
  345.         // used in interpolation
  346.         // this imply the ancestor of min/max at (n, m) is one of
  347.         // (2n, m), (2n+1, m), (2n+2, m), (2n, m+1), (2n+1, m+1), (2n+2, m+1)
  348.         int selectedI = levelI;
  349.         int selectedJ = 2 * levelJ;
  350.         double selectedElevation = Double.NaN;
  351.         for (int n = 2 * levelJ; n < 2 * levelJ + 3; ++n) {
  352.             if (n < getLongitudeColumns()) {
  353.                 for (int m = levelI; m < levelI + 2; ++m) {
  354.                     if (m < getLatitudeRows()) {
  355.                         final double elevation = raw[m * getLongitudeColumns() + n];
  356.                         if (selector.selectFirst(elevation, selectedElevation)) {
  357.                             selectedI         = m;
  358.                             selectedJ         = n;
  359.                             selectedElevation = elevation;
  360.                         }
  361.                     }
  362.                 }
  363.             }
  364.         }

  366.         return new int[] {
  367.             selectedI, selectedJ
  368.         };

  370.     }

  372.     /** Get the deepest level at which two cells are merged in the same min/max sub-tile.
  373.      * @param i1 row index of first cell
  374.      * @param j1 column index of first cell
  375.      * @param i2 row index of second cell
  376.      * @param j2 column index of second cell
  377.      * @return deepest level at which two cells are merged in the same min/max sub-tile,
  378.      * or -1 if they are never merged in the same sub-tile
  379.      * @see #getLevels()
  380.      * @see #getMinElevation(int, int, int)
  381.      * @see #getMaxElevation(int, int, int)
  382.      */
  383.     public int getMergeLevel(final int i1, final int j1, final int i2, final int j2) {

  385.         int largest = -1;

  387.         for (int level = 0; level < start.length; ++level) {
  388.             // compute indices in level merged array
  389.             final int k        = start.length - level;
  390.             final int rowShift = k / 2;
  391.             final int colShift = (k + 1) / 2;
  392.             final int levelI1  = i1 >> rowShift;
  393.             final int levelJ1  = j1 >> colShift;
  394.             final int levelI2  = i2 >> rowShift;
  395.             final int levelJ2  = j2 >> colShift;
  396.             if (levelI1 != levelI2 || levelJ1 != levelJ2) {
  397.                 return largest;
  398.             }
  399.             largest = level;
  400.         }

  402.         return largest;

  404.     }

  406.     /** Get the index of sub-tiles start rows crossed.
  407.      * <p>
  408.      * When going from one row to another row at some tree level,
  409.      * we cross sub-tiles boundaries. This method returns the index
  410.      * of these boundaries.
  411.      * </p>
  412.      * @param row1 starting row
  413.      * @param row2 ending row
  414.      * @param level tree level
  415.      * @return indices of rows crossed at sub-tiles boundaries, in crossing order,
  416.      * the endpoints <em>are</em> included (i.e. if {@code row1} or {@code row2} are
  417.      * boundary rows, they will be in returned array)
  418.      */
  419.     public int[] getCrossedBoundaryRows(final int row1, final int row2, final int level) {

  421.         // number of rows in each sub-tile
  422.         final int rows = 1 << ((start.length - level) / 2);

  424.         // build the crossings in ascending order
  425.         final int min = FastMath.min(row1, row2);
  426.         final int max = FastMath.max(row1, row2) + 1;
  427.         return buildCrossings((min + rows - 1) - ((min + rows - 1) % rows), max, rows,
  428.                               row1 <= row2);

  430.     }

  432.     /** Get the index of sub-tiles start columns crossed.
  433.      * <p>
  434.      * When going from one column to another column at some tree level,
  435.      * we cross sub-tiles boundaries. This method returns the index
  436.      * of these boundaries.
  437.      * </p>
  438.      * @param column1 starting column
  439.      * @param column2 ending column (excluded)
  440.      * @param level tree level
  441.      * @return indices of columns crossed at sub-tiles boundaries, in crossing order,
  442.      * the endpoints <em>are</em> included (i.e. if {@code column1} or {@code column2} are
  443.      * boundary columns, they will be in returned array)
  444.      */
  445.     public int[] getCrossedBoundaryColumns(final int column1, final int column2, final int level) {

  447.         // number of columns in each sub-tile
  448.         final int columns  = 1 << ((start.length + 1 - level) / 2);

  450.         // build the crossings in ascending order
  451.         final int min = FastMath.min(column1, column2);
  452.         final int max = FastMath.max(column1, column2) + 1;
  453.         return buildCrossings((min + columns - 1) - ((min + columns - 1) % columns), max, columns,
  454.                               column1 <= column2);

  456.     }

  458.     /** Build crossings arrays.
  459.      * @param begin begin crossing index
  460.      * @param end end crossing index (excluded, if equal to begin, the array is empty)
  461.      * @param step crossing step
  462.      * @param ascending if true, the crossings must be in ascending order
  463.      * @return indices of rows or columns crossed at sub-tiles boundaries, in crossing order
  464.      */
  465.     private int[] buildCrossings(final int begin, final int end, final int step, final boolean ascending) {

  467.         // allocate array
  468.         final int n = FastMath.max(0, (end - begin + step + ((step > 0) ? -1 : +1)) / step);
  469.         final int[] crossings = new int[n];

  471.         // fill it up
  472.         int crossing = begin;
  473.         if (ascending) {
  474.             for (int i = 0; i < crossings.length; ++i) {
  475.                 crossings[i] = crossing;
  476.                 crossing += step;
  477.             }
  478.         } else {
  479.             for (int i = 0; i < crossings.length; ++i) {
  480.                 crossings[crossings.length - 1 - i] = crossing;
  481.                 crossing += step;
  482.             }
  483.         }

  485.         return crossings;

  487.     }

  489.     /** Check if the merging operation between level and level-1 is a column merging.
  490.      * @param level level to check
  491.      * @return true if the merging operation between level and level-1 is a column
  492.      * merging, false if is a row merging
  493.      */
  494.     public boolean isColumnMerging(final int level) {
  495.         return (level & 0x1) == (start.length & 0x1);
  496.     }

  498.     /** Recursive setting of tree levels.
  499.      * <p>
  500.      * The following algorithms works for any array shape, even with
  501.      * rows or columns which are not powers of 2 or with one
  502.      * dimension much larger than the other. As an example, starting
  503.      * from a 107 ⨉ 19 array, we get the following 9 levels, for a
  504.      * total of 2187 elements in each tree:
  505.      * </p>
  506.      * <p>
  507.      * <table border="0">
  508.      * <tr BGCOLOR="#EEEEFF">
  509.      *     <td>Level</td>   <td>Dimension</td>  <td>Start index</td>  <td>End index (inclusive)</td></tr>
  510.      * <tr>   <td>0</td>     <td>  7 ⨉  1</td>       <td>   0</td>        <td>  6</td> </tr>
  511.      * <tr>   <td>1</td>     <td>  7 ⨉  2</td>       <td>   7</td>        <td> 20</td> </tr>
  512.      * <tr>   <td>2</td>     <td> 14 ⨉  2</td>       <td>  21</td>        <td> 48</td> </tr>
  513.      * <tr>   <td>3</td>     <td> 14 ⨉  3</td>       <td>  49</td>        <td> 90</td> </tr>
  514.      * <tr>   <td>4</td>     <td> 27 ⨉  3</td>       <td>  91</td>        <td>171</td> </tr>
  515.      * <tr>   <td>5</td>     <td> 27 ⨉  5</td>       <td> 172</td>        <td>306</td> </tr>
  516.      * <tr>   <td>6</td>     <td> 54 ⨉  5</td>       <td> 307</td>        <td>576</td> </tr>
  517.      * <tr>   <td>7</td>     <td> 54 ⨉ 10</td>      <td> 577</td>        <td>1116</td> </tr>
  518.      * <tr>   <td>8</td>     <td>107 ⨉ 10</td>      <td>1117</td>        <td>2186</td> </tr>
  519.      * </table>
  520.      * </p>
  521.      * @param stage number of merging stages
  522.      * @param stageRows number of rows at current stage
  523.      * @param stageColumns number of columns at current stage
  524.      * @return size cumulative size from root to current level
  525.      */
  526.     private int setLevels(final int stage, final int stageRows, final int stageColumns) {

  528.         if (stageRows == 1 || stageColumns == 1) {
  529.             // we have found root, stop recursion
  530.             start   = new int[stage];
  531.             if (stage > 0) {
  532.                 start[0]   = 0;
  533.             }
  534.             return stageRows * stageColumns;
  535.         }

  537.         final int size;
  538.         if ((stage & 0x1) == 0) {
  539.             // columns merging
  540.             size = setLevels(stage + 1, stageRows, (stageColumns + 1) / 2);
  541.         } else {
  542.             // rows merging
  543.             size = setLevels(stage + 1, (stageRows + 1) / 2, stageColumns);
  544.         }

  546.         if (stage > 0) {
  547.             // store current level characteristics
  548.             start[start.length     - stage] = size;
  549.             return size + stageRows * stageColumns;
  550.         } else {
  551.             // we don't count the elements at stage 0 as they are not stored in the
  552.             // min/max trees (they correspond to the raw elevation, without merging)
  553.             return size;
  554.         }

  556.     }

  558.     /** Preprocess recursive application of a function.
  559.      * <p>
  560.      * At start, the min/max should be computed for each cell using the four corners values.
  561.      * </p>
  562.      * @param preprocessed preprocessed array to fill up
  563.      * @param elevations raw elevations te preprocess
  564.      * @param nbRows number of rows
  565.      * @param nbCols number of columns
  566.      * @param selector selector to use
  567.      */
  568.     private void preprocess(final double[] preprocessed, final double[] elevations,
  569.                             final int nbRows, final int nbCols,
  570.                             final Selector selector) {

  572.         int k = 0;

  574.         for (int i = 0; i < nbRows - 1; ++i) {

  576.             // regular elements with both a column at right and a row below
  577.             for (int j = 0; j < nbCols - 1; ++j) {
  578.                 preprocessed[k] = selector.select(selector.select(elevations[k],          elevations[k + 1]),
  579.                                                   selector.select(elevations[k + nbCols], elevations[k + nbCols + 1]));
  580.                 k++;
  581.             }

  583.             // last column elements, lacking a right column
  584.             preprocessed[k] = selector.select(elevations[k], elevations[k + nbCols]);
  585.             k++;

  587.         }

  589.         // last row elements, lacking a below row
  590.         for (int j = 0; j < nbCols - 1; ++j) {
  591.             preprocessed[k] = selector.select(elevations[k], elevations[k + 1]);
  592.             k++;
  593.         }

  595.         // last element
  596.         preprocessed[k] = elevations[k];

  598.     }

  600.     /** Recursive application of a function.
  601.      * @param tree tree to fill-up with the recursive applications
  602.      * @param level current level
  603.      * @param levelRows number of rows at current level
  604.      * @param levelColumns number of columns at current level
  605.      * @param selector to apply
  606.      * @param base base array from which function arguments are drawn
  607.      * @param first index of the first element to consider in base array
  608.      */
  609.     private void applyRecursively(final double[] tree,
  610.                                   final int level, final int levelRows, final int levelColumns,
  611.                                   final Selector selector,
  612.                                   final double[] base, final int first) {

  614.         if (isColumnMerging(level + 1)) {

  616.             // merge columns pairs
  617.             int           iTree       = start[level];
  618.             int           iBase       = first;
  619.             final int     nextColumns = (levelColumns + 1) / 2;
  620.             final boolean odd         = (levelColumns & 0x1) != 0;
  621.             final int     jEnd        = odd ? nextColumns - 1 : nextColumns;
  622.             for (int i = 0; i < levelRows; ++i) {

  624.                 // regular pairs
  625.                 for (int j = 0; j < jEnd; ++j) {
  626.                     tree[iTree++] = selector.select(base[iBase], base[iBase + 1]);
  627.                     iBase += 2;
  628.                 }

  630.                 if (odd) {
  631.                     // last column
  632.                     tree[iTree++] = base[iBase++];
  633.                 }


  636.             }

  638.             if (level > 0) {
  639.                 applyRecursively(tree, level - 1, levelRows, nextColumns, selector, tree, start[level]);
  640.             }

  642.         } else {

  644.             // merge rows pairs
  645.             int           iTree    = start[level];
  646.             int           iBase    = first;
  647.             final int     nextRows = (levelRows + 1) / 2;
  648.             final boolean odd      = (levelRows & 0x1) != 0;
  649.             final int     iEnd     = odd ? nextRows - 1 : nextRows;

  651.             // regular pairs
  652.             for (int i = 0; i < iEnd; ++i) {

  654.                 for (int j = 0; j < levelColumns; ++j) {
  655.                     tree[iTree++] = selector.select(base[iBase], base[iBase + levelColumns]);
  656.                     iBase++;
  657.                 }
  658.                 iBase += levelColumns;

  660.             }

  662.             if (odd) {
  663.                 // last row
  664.                 System.arraycopy(base, iBase, tree, iTree, levelColumns);
  665.             }

  667.             if (level > 0) {
  668.                 applyRecursively(tree, level - 1, nextRows, levelColumns, selector, tree, start[level]);
  669.             }

  671.         }
  672.     }

  674. }
Created with JaCoCo 0.8.7.202105040129