1 /* Copyright 2013-2019 CS Systèmes d'Information
2 * Licensed to CS Systèmes d'Information (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;
18
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;
25
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 {
75
76 /** Raw elevations. */
77 private double[] raw;
78
79 /** Min kd-tree. */
80 private double[] minTree;
81
82 /** Max kd-tree. */
83 private double[] maxTree;
84
85 /** Start indices of tree levels. */
86 private int[] start;
87
88 /** Simple constructor.
89 * <p>
90 * Creates an empty tile.
91 * </p>
92 */
93 MinMaxTreeTile() {
94 }
95
96 /** {@inheritDoc} */
97 @Override
98 protected void processUpdatedElevation(final double[] elevations) {
99
100 raw = elevations;
101
102 final int nbRows = getLatitudeRows();
103 final int nbCols = getLongitudeColumns();
104
105 // set up the levels
106 final int size = setLevels(0, nbRows, nbCols);
107 minTree = new double[size];
108 maxTree = new double[size];
109
110 // compute min/max trees
111 if (start.length > 0) {
112
113 final double[] preprocessed = new double[raw.length];
114
115 preprocess(preprocessed, raw, nbRows, nbCols, MinSelector.getInstance());
116 applyRecursively(minTree, start.length - 1, nbRows, nbCols, MinSelector.getInstance(), preprocessed, 0);
117
118 preprocess(preprocessed, raw, nbRows, nbCols, MaxSelector.getInstance());
119 applyRecursively(maxTree, start.length - 1, nbRows, nbCols, MaxSelector.getInstance(), preprocessed, 0);
120
121 }
122
123 }
124
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 }
134
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) {
167
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);
175
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 }
190
191 return minTree[start[level] + levelI * levelC + levelJ];
192
193 }
194
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) {
227
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);
235
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 }
250
251 return maxTree[start[level] + levelI * levelC + levelJ];
252
253 }
254
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 }
271
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 }
288
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) {
299
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);
307
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) {
310
311 if (isColumnMerging(l)) {
312
313 --colShift;
314 levelC = 1 + ((getLongitudeColumns() - 1) >> colShift);
315 levelJ = levelJ << 1;
316
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 }
324
325 } else {
326
327 --rowShift;
328 levelR = 1 + ((getLatitudeRows() - 1) >> rowShift);
329 levelI = levelI << 1;
330
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 }
338
339 }
340
341 }
342
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 }
365
366 return new int[] {
367 selectedI, selectedJ
368 };
369
370 }
371
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) {
384
385 int largest = -1;
386
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 }
401
402 return largest;
403
404 }
405
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) {
420
421 // number of rows in each sub-tile
422 final int rows = 1 << ((start.length - level) / 2);
423
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);
429
430 }
431
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) {
446
447 // number of columns in each sub-tile
448 final int columns = 1 << ((start.length + 1 - level) / 2);
449
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);
455
456 }
457
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) {
466
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];
470
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 }
484
485 return crossings;
486
487 }
488
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 }
497
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) {
527
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 }
536
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 }
545
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 }
555
556 }
557
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) {
571
572 int k = 0;
573
574 for (int i = 0; i < nbRows - 1; ++i) {
575
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 }
582
583 // last column elements, lacking a right column
584 preprocessed[k] = selector.select(elevations[k], elevations[k + nbCols]);
585 k++;
586
587 }
588
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 }
594
595 // last element
596 preprocessed[k] = elevations[k];
597
598 }
599
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) {
613
614 if (isColumnMerging(level + 1)) {
615
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) {
623
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 }
629
630 if (odd) {
631 // last column
632 tree[iTree++] = base[iBase++];
633 }
634
635
636 }
637
638 if (level > 0) {
639 applyRecursively(tree, level - 1, levelRows, nextColumns, selector, tree, start[level]);
640 }
641
642 } else {
643
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;
650
651 // regular pairs
652 for (int i = 0; i < iEnd; ++i) {
653
654 for (int j = 0; j < levelColumns; ++j) {
655 tree[iTree++] = selector.select(base[iBase], base[iBase + levelColumns]);
656 iBase++;
657 }
658 iBase += levelColumns;
659
660 }
661
662 if (odd) {
663 // last row
664 System.arraycopy(base, iBase, tree, iTree, levelColumns);
665 }
666
667 if (level > 0) {
668 applyRecursively(tree, level - 1, nextRows, levelColumns, selector, tree, start[level]);
669 }
670
671 }
672 }
673
674 }