1 /* Copyright 2013-2025 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;
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>
58 * <caption>"min/max kd-tree levels for a 107 x19 raw tile "</caption>
59 * <thead>
60 * <tr>
61 * <th>Level</th>
62 * <th>Number of sub-tiles</th>
63 * <th>Regular sub-tiles dimension</th>
64 * </tr>
65 * </thead>
66 * <tbody>
67 * <tr>
68 * <td>8</td>
69 * <td>107 ⨉ 10</td>
70 * <td>1 ⨉ 2</td>
71 * </tr>
72 * <tr>
73 * <td>7</td>
74 * <td>54 ⨉ 10</td>
75 * <td>2 ⨉ 2</td>
76 * </tr>
77 * <tr>
78 * <td>6</td>
79 * <td>54 ⨉ 5</td>
80 * <td>2 ⨉ 4</td>
81 * </tr>
82 * <tr>
83 * <td>5</td>
84 * <td>27 ⨉ 5</td>
85 * <td>4 ⨉ 4</td>
86 * </tr>
87 * <tr>
88 * <td>4</td>
89 * <td>27 ⨉ 3</td>
90 * <td>4 ⨉ 8</td>
91 * </tr>
92 * <tr>
93 * <td>3</td>
94 * <td>14 ⨉ 3</td>
95 * <td>8 ⨉ 8</td>
96 * </tr>
97 * <tr>
98 * <td>2</td>
99 * <td>14 ⨉ 2</td>
100 * <td>8 ⨉ 16</td>
101 * </tr>
102 * <tr>
103 * <td>1</td>
104 * <td>7 ⨉ 2</td>
105 * <td>16 ⨉ 16</td>
106 * </tr>
107 * <tr>
108 * <td>0</td>
109 * <td>7 ⨉ 1</td>
110 * <td>16 ⨉ 32</td>
111 * </tr>
112 * </tbody>
113 *</table>
114 *
115 * @see MinMaxTreeTileFactory
116 * @author Luc Maisonobe
117 */
118 public class MinMaxTreeTile extends SimpleTile {
119
120 /** Raw elevations. */
121 private double[] raw;
122
123 /** Min kd-tree. */
124 private double[] minTree;
125
126 /** Max kd-tree. */
127 private double[] maxTree;
128
129 /** Start indices of tree levels. */
130 private int[] start;
131
132 /** Simple constructor.
133 * <p>
134 * Creates an empty tile.
135 * </p>
136 */
137 protected MinMaxTreeTile() {
138 }
139
140 /** {@inheritDoc} */
141 @Override
142 protected void processUpdatedElevation(final double[] elevations) {
143
144 raw = elevations;
145
146 final int nbRows = getLatitudeRows();
147 final int nbCols = getLongitudeColumns();
148
149 // set up the levels
150 final int size = setLevels(0, nbRows, nbCols);
151 minTree = new double[size];
152 maxTree = new double[size];
153
154 // compute min/max trees
155 if (start.length > 0) {
156
157 final double[] preprocessed = new double[raw.length];
158
159 preprocess(preprocessed, raw, nbRows, nbCols, MinSelector.getInstance());
160 applyRecursively(minTree, start.length - 1, nbRows, nbCols, MinSelector.getInstance(), preprocessed, 0);
161
162 preprocess(preprocessed, raw, nbRows, nbCols, MaxSelector.getInstance());
163 applyRecursively(maxTree, start.length - 1, nbRows, nbCols, MaxSelector.getInstance(), preprocessed, 0);
164
165 }
166
167 }
168
169 /** Get the number of kd-tree levels (not counting raw elevations).
170 * @return number of kd-tree levels
171 * @see #getMinElevation(int, int, int)
172 * @see #getMaxElevation(int, int, int)
173 * @see #getMergeLevel(int, int, int, int)
174 */
175 public int getLevels() {
176 return start.length;
177 }
178
179 /** Get the minimum elevation at some level tree.
180 * <p>
181 * Note that the min elevation is <em>not</em> computed
182 * only at cell center, but considering that it is interpolated
183 * considering also Eastwards and Northwards neighbors, and extends
184 * up to the center of these neighbors. As an example, lets consider
185 * four neighboring cells in some Digital Elevation Model:
186 *
187 * <table>
188 * <caption>"four neighboring cells"</caption>
189 *<thead>
190 * <tr><th>j/i</th> <th>i</th> <th>i+1</th></tr>
191 *</thead>
192 *<tbody>
193 * <tr> <th>j+1</th><td>11</td><td>12</td> </tr>
194 * <tr> <th>j</th> <td>10</td> <td>11</td> </tr>
195 * </tbody>
196 *</table>
197 * When we interpolate elevation at a point located slightly South-West
198 * to the center of the (i+1, j+1) cell, we use all four cells in the
199 * interpolation, and we will get a result very close to 10 if we start
200 * close to (i+1, j+1) cell center. As the min value for this interpolation
201 * is stored at (i, j) indices, this implies that {@code getMinElevation(i,
202 * j, l)} must return 10 if l is chosen such that the sub-tile at
203 * tree level l includes cell (i,j) but not cell (i+1, j+1). In other words,
204 * interpolation implies sub-tile boundaries are overshoot by one column to
205 * the East and one row to the North when computing min.
206 *
207 * @param i row index of the cell
208 * @param j column index of the cell
209 * @param level tree level
210 * @return minimum value that can be reached when interpolating elevation
211 * in the sub-tile
212 * @see #getLevels()
213 * @see #getMaxElevation(int, int, int)
214 * @see #getMergeLevel(int, int, int, int)
215 */
216 public double getMinElevation(final int i, final int j, final int level) {
217
218 // compute indices in level merged array
219 final int k = start.length - level;
220 final int rowShift = k / 2;
221 final int colShift = (k + 1) / 2;
222 final int levelI = i >> rowShift;
223 final int levelJ = j >> colShift;
224 final int levelC = 1 + ((getLongitudeColumns() - 1) >> colShift);
225
226 if (DumpManager.isActive()) {
227 final int[] min = locateMin(i, j, level);
228 final int index = min[0] * getLongitudeColumns() + min[1];
229 DumpManager.dumpTileCell(this, min[0], min[1], raw[index]);
230 if (index + getLongitudeColumns() < raw.length) {
231 DumpManager.dumpTileCell(this, min[0] + 1, min[1], raw[index + getLongitudeColumns()]);
232 }
233 if (index + 1 < raw.length) {
234 DumpManager.dumpTileCell(this, min[0], min[1] + 1, raw[index + 1]);
235 }
236 if (index + getLongitudeColumns() + 1 < raw.length) {
237 DumpManager.dumpTileCell(this, min[0] + 1, min[1] + 1, raw[index + getLongitudeColumns() + 1]);
238 }
239 }
240
241 return minTree[start[level] + levelI * levelC + levelJ];
242
243 }
244
245 /** Get the maximum elevation at some level tree.
246 * <p>
247 * Note that the max elevation is <em>not</em> computed
248 * only at cell center, but considering that it is interpolated
249 * considering also Eastwards and Northwards neighbors, and extends
250 * up to the center of these neighbors. As an example, lets consider
251 * four neighboring cells in some Digital Elevation Model:
252 * <table>
253 * <caption>"four neighboring cells"</caption>
254 *<thead>
255 * <tr><th>j/i</th> <th>i</th> <th>i+1</th></tr>
256 *</thead>
257 *<tbody>
258 * <tr> <th>j+1</th><td>11</td><td>10</td> </tr>
259 * <tr> <th>j</th> <td>12</td> <td>11</td> </tr>
260 * </tbody>
261 *</table>
262 * When we interpolate elevation at a point located slightly South-West
263 * to the center of the (i+1, j+1) cell, we use all four cells in the
264 * interpolation, and we will get a result very close to 12 if we start
265 * close to (i+1, j+1) cell center. As the max value for this interpolation
266 * is stored at (i, j) indices, this implies that {@code getMaxElevation(i,
267 * j, l)} must return 12 if l is chosen such that the sub-tile at
268 * tree level l includes cell (i,j) but not cell (i+1, j+1). In other words,
269 * interpolation implies sub-tile boundaries are overshoot by one column to
270 * the East and one row to the North when computing max.
271 *
272 * @param i row index of the cell
273 * @param j column index of the cell
274 * @param level tree level
275 * @return maximum value that can be reached when interpolating elevation
276 * in the sub-tile
277 * @see #getLevels()
278 * @see #getMinElevation(int, int, int)
279 * @see #getMergeLevel(int, int, int, int)
280 */
281 public double getMaxElevation(final int i, final int j, final int level) {
282
283 // compute indices in level merged array
284 final int k = start.length - level;
285 final int rowShift = k / 2;
286 final int colShift = (k + 1) / 2;
287 final int levelI = i >> rowShift;
288 final int levelJ = j >> colShift;
289 final int levelC = 1 + ((getLongitudeColumns() - 1) >> colShift);
290
291 if (DumpManager.isActive()) {
292 final int[] max = locateMax(i, j, level);
293 final int index = max[0] * getLongitudeColumns() + max[1];
294 DumpManager.dumpTileCell(this, max[0], max[1], raw[index]);
295 if (index + getLongitudeColumns() < raw.length) {
296 DumpManager.dumpTileCell(this, max[0] + 1, max[1], raw[index + getLongitudeColumns()]);
297 }
298 if (index + 1 < raw.length) {
299 DumpManager.dumpTileCell(this, max[0], max[1] + 1, raw[index + 1]);
300 }
301 if (index + getLongitudeColumns() + 1 < raw.length) {
302 DumpManager.dumpTileCell(this, max[0] + 1, max[1] + 1, raw[index + getLongitudeColumns() + 1]);
303 }
304 }
305
306 return maxTree[start[level] + levelI * levelC + levelJ];
307
308 }
309
310 /** Locate the cell at which min elevation is reached for a specified level.
311 * <p>
312 * Min is computed with respect to the continuous interpolated elevation, which
313 * takes four neighboring cells into account. This implies that the cell at which
314 * min value is reached for some level is either within the sub-tile for this level,
315 * or in some case it may be one column outside to the East or one row outside to
316 * the North. See {@link #getMinElevation()} for a more complete explanation.
317 * </p>
318 * @param i row index of the cell
319 * @param j column index of the cell
320 * @param level tree level of the sub-tile considered
321 * @return row/column indices of the cell at which min elevation is reached
322 */
323 public int[] locateMin(final int i, final int j, final int level) {
324 return locateMinMax(i, j, level, MinSelector.getInstance(), minTree);
325 }
326
327 /** Locate the cell at which max elevation is reached for a specified level.
328 * <p>
329 * Max is computed with respect to the continuous interpolated elevation, which
330 * takes four neighboring cells into account. This implies that the cell at which
331 * max value is reached for some level is either within the sub-tile for this level,
332 * or in some case it may be one column outside to the East or one row outside to
333 * the North. See {@link #getMaxElevation()} for a more complete explanation.
334 * </p>
335 * @param i row index of the cell
336 * @param j column index of the cell
337 * @param level tree level of the sub-tile considered
338 * @return row/column indices of the cell at which min elevation is reached
339 */
340 public int[] locateMax(final int i, final int j, final int level) {
341 return locateMinMax(i, j, level, MaxSelector.getInstance(), maxTree);
342 }
343
344 /** Locate the cell at which min/max elevation is reached for a specified level.
345 * @param i row index of the cell
346 * @param j column index of the cell
347 * @param level tree level of the sub-tile considered
348 * @param selector min/max selector to use
349 * @param tree min/max tree to use
350 * @return row/column indices of the cell at which min/max elevation is reached
351 */
352 private int[] locateMinMax(final int i, final int j, final int level,
353 final Selector selector, final double[] tree) {
354
355 final int k = start.length - level;
356 int rowShift = k / 2;
357 int colShift = (k + 1) / 2;
358 int levelI = i >> rowShift;
359 int levelJ = j >> colShift;
360 int levelR = 1 + ((getLatitudeRows() - 1) >> rowShift);
361 int levelC = 1 + ((getLongitudeColumns() - 1) >> colShift);
362
363 // track the cell ancestors from merged tree at specified level up to tree at level 1
364 for (int l = level + 1; l < start.length; ++l) {
365
366 if (isColumnMerging(l)) {
367
368 --colShift;
369 levelC = 1 + ((getLongitudeColumns() - 1) >> colShift);
370 levelJ = levelJ << 1;
371
372 if (levelJ + 1 < levelC) {
373 // the cell results from a regular merging of two columns
374 if (selector.selectFirst(tree[start[l] + levelI * levelC + levelJ + 1],
375 tree[start[l] + levelI * levelC + levelJ])) {
376 levelJ++;
377 }
378 }
379
380 } else {
381
382 --rowShift;
383 levelR = 1 + ((getLatitudeRows() - 1) >> rowShift);
384 levelI = levelI << 1;
385
386 if (levelI + 1 < levelR) {
387 // the cell results from a regular merging of two rows
388 if (selector.selectFirst(tree[start[l] + (levelI + 1) * levelC + levelJ],
389 tree[start[l] + levelI * levelC + levelJ])) {
390 levelI++;
391 }
392 }
393
394 }
395
396 }
397
398 // we are now at first merge level, which always results from a column merge
399 // or pre-processed data, which themselves result from merging four cells
400 // used in interpolation
401 // this imply the ancestor of min/max at (n, m) is one of
402 // (2n, m), (2n+1, m), (2n+2, m), (2n, m+1), (2n+1, m+1), (2n+2, m+1)
403 int selectedI = levelI;
404 int selectedJ = 2 * levelJ;
405 double selectedElevation = Double.NaN;
406 for (int n = 2 * levelJ; n < 2 * levelJ + 3; ++n) {
407 if (n < getLongitudeColumns()) {
408 for (int m = levelI; m < levelI + 2; ++m) {
409 if (m < getLatitudeRows()) {
410 final double elevation = raw[m * getLongitudeColumns() + n];
411 if (selector.selectFirst(elevation, selectedElevation)) {
412 selectedI = m;
413 selectedJ = n;
414 selectedElevation = elevation;
415 }
416 }
417 }
418 }
419 }
420
421 return new int[] {
422 selectedI, selectedJ
423 };
424
425 }
426
427 /** Get the deepest level at which two cells are merged in the same min/max sub-tile.
428 * @param i1 row index of first cell
429 * @param j1 column index of first cell
430 * @param i2 row index of second cell
431 * @param j2 column index of second cell
432 * @return deepest level at which two cells are merged in the same min/max sub-tile,
433 * or -1 if they are never merged in the same sub-tile
434 * @see #getLevels()
435 * @see #getMinElevation(int, int, int)
436 * @see #getMaxElevation(int, int, int)
437 */
438 public int getMergeLevel(final int i1, final int j1, final int i2, final int j2) {
439
440 int largest = -1;
441
442 for (int level = 0; level < start.length; ++level) {
443 // compute indices in level merged array
444 final int k = start.length - level;
445 final int rowShift = k / 2;
446 final int colShift = (k + 1) / 2;
447 final int levelI1 = i1 >> rowShift;
448 final int levelJ1 = j1 >> colShift;
449 final int levelI2 = i2 >> rowShift;
450 final int levelJ2 = j2 >> colShift;
451 if (levelI1 != levelI2 || levelJ1 != levelJ2) {
452 return largest;
453 }
454 largest = level;
455 }
456
457 return largest;
458
459 }
460
461 /** Get the index of sub-tiles start rows crossed.
462 * <p>
463 * When going from one row to another row at some tree level,
464 * we cross sub-tiles boundaries. This method returns the index
465 * of these boundaries.
466 * </p>
467 * @param row1 starting row
468 * @param row2 ending row
469 * @param level tree level
470 * @return indices of rows crossed at sub-tiles boundaries, in crossing order,
471 * the endpoints <em>are</em> included (i.e. if {@code row1} or {@code row2} are
472 * boundary rows, they will be in returned array)
473 */
474 public int[] getCrossedBoundaryRows(final int row1, final int row2, final int level) {
475
476 // number of rows in each sub-tile
477 final int rows = 1 << ((start.length - level) / 2);
478
479 // build the crossings in ascending order
480 final int min = FastMath.min(row1, row2);
481 final int max = FastMath.max(row1, row2) + 1;
482 return buildCrossings((min + rows - 1) - ((min + rows - 1) % rows), max, rows,
483 row1 <= row2);
484
485 }
486
487 /** Get the index of sub-tiles start columns crossed.
488 * <p>
489 * When going from one column to another column at some tree level,
490 * we cross sub-tiles boundaries. This method returns the index
491 * of these boundaries.
492 * </p>
493 * @param column1 starting column
494 * @param column2 ending column (excluded)
495 * @param level tree level
496 * @return indices of columns crossed at sub-tiles boundaries, in crossing order,
497 * the endpoints <em>are</em> included (i.e. if {@code column1} or {@code column2} are
498 * boundary columns, they will be in returned array)
499 */
500 public int[] getCrossedBoundaryColumns(final int column1, final int column2, final int level) {
501
502 // number of columns in each sub-tile
503 final int columns = 1 << ((start.length + 1 - level) / 2);
504
505 // build the crossings in ascending order
506 final int min = FastMath.min(column1, column2);
507 final int max = FastMath.max(column1, column2) + 1;
508 return buildCrossings((min + columns - 1) - ((min + columns - 1) % columns), max, columns,
509 column1 <= column2);
510
511 }
512
513 /** Build crossings arrays.
514 * @param begin begin crossing index
515 * @param end end crossing index (excluded, if equal to begin, the array is empty)
516 * @param step crossing step
517 * @param ascending if true, the crossings must be in ascending order
518 * @return indices of rows or columns crossed at sub-tiles boundaries, in crossing order
519 */
520 private int[] buildCrossings(final int begin, final int end, final int step, final boolean ascending) {
521
522 // allocate array
523 final int n = FastMath.max(0, (end - begin + step + ((step > 0) ? -1 : +1)) / step);
524 final int[] crossings = new int[n];
525
526 // fill it up
527 int crossing = begin;
528 if (ascending) {
529 for (int i = 0; i < crossings.length; ++i) {
530 crossings[i] = crossing;
531 crossing += step;
532 }
533 } else {
534 for (int i = 0; i < crossings.length; ++i) {
535 crossings[crossings.length - 1 - i] = crossing;
536 crossing += step;
537 }
538 }
539
540 return crossings;
541
542 }
543
544 /** Check if the merging operation between level and level-1 is a column merging.
545 * @param level level to check
546 * @return true if the merging operation between level and level-1 is a column
547 * merging, false if is a row merging
548 */
549 public boolean isColumnMerging(final int level) {
550 return (level & 0x1) == (start.length & 0x1);
551 }
552
553 /** Recursive setting of tree levels.
554 * <p>
555 * The following algorithms works for any array shape, even with
556 * rows or columns which are not powers of 2 or with one
557 * dimension much larger than the other. As an example, starting
558 * from a 107 ⨉ 19 array, we get the following 9 levels, for a
559 * total of 2187 elements in each tree:
560 * </p>
561 * <p>
562 * <table border="0">
563 * <tr BGCOLOR="#EEEEFF">
564 * <td>Level</td> <td>Dimension</td> <td>Start index</td> <td>End index (inclusive)</td></tr>
565 * <tr> <td>0</td> <td> 7 ⨉ 1</td> <td> 0</td> <td> 6</td> </tr>
566 * <tr> <td>1</td> <td> 7 ⨉ 2</td> <td> 7</td> <td> 20</td> </tr>
567 * <tr> <td>2</td> <td> 14 ⨉ 2</td> <td> 21</td> <td> 48</td> </tr>
568 * <tr> <td>3</td> <td> 14 ⨉ 3</td> <td> 49</td> <td> 90</td> </tr>
569 * <tr> <td>4</td> <td> 27 ⨉ 3</td> <td> 91</td> <td>171</td> </tr>
570 * <tr> <td>5</td> <td> 27 ⨉ 5</td> <td> 172</td> <td>306</td> </tr>
571 * <tr> <td>6</td> <td> 54 ⨉ 5</td> <td> 307</td> <td>576</td> </tr>
572 * <tr> <td>7</td> <td> 54 ⨉ 10</td> <td> 577</td> <td>1116</td> </tr>
573 * <tr> <td>8</td> <td>107 ⨉ 10</td> <td>1117</td> <td>2186</td> </tr>
574 * </table>
575 * </p>
576 * @param stage number of merging stages
577 * @param stageRows number of rows at current stage
578 * @param stageColumns number of columns at current stage
579 * @return size cumulative size from root to current level
580 */
581 private int setLevels(final int stage, final int stageRows, final int stageColumns) {
582
583 if (stageRows == 1 || stageColumns == 1) {
584 // we have found root, stop recursion
585 start = new int[stage];
586 if (stage > 0) {
587 start[0] = 0;
588 }
589 return stageRows * stageColumns;
590 }
591
592 final int size;
593 if ((stage & 0x1) == 0) {
594 // columns merging
595 size = setLevels(stage + 1, stageRows, (stageColumns + 1) / 2);
596 } else {
597 // rows merging
598 size = setLevels(stage + 1, (stageRows + 1) / 2, stageColumns);
599 }
600
601 if (stage > 0) {
602 // store current level characteristics
603 start[start.length - stage] = size;
604 return size + stageRows * stageColumns;
605 } else {
606 // we don't count the elements at stage 0 as they are not stored in the
607 // min/max trees (they correspond to the raw elevation, without merging)
608 return size;
609 }
610
611 }
612
613 /** Preprocess recursive application of a function.
614 * <p>
615 * At start, the min/max should be computed for each cell using the four corners values.
616 * </p>
617 * @param preprocessed preprocessed array to fill up
618 * @param elevations raw elevations te preprocess
619 * @param nbRows number of rows
620 * @param nbCols number of columns
621 * @param selector selector to use
622 */
623 private void preprocess(final double[] preprocessed, final double[] elevations,
624 final int nbRows, final int nbCols,
625 final Selector selector) {
626
627 int k = 0;
628
629 for (int i = 0; i < nbRows - 1; ++i) {
630
631 // regular elements with both a column at right and a row below
632 for (int j = 0; j < nbCols - 1; ++j) {
633 preprocessed[k] = selector.select(selector.select(elevations[k], elevations[k + 1]),
634 selector.select(elevations[k + nbCols], elevations[k + nbCols + 1]));
635 k++;
636 }
637
638 // last column elements, lacking a right column
639 preprocessed[k] = selector.select(elevations[k], elevations[k + nbCols]);
640 k++;
641
642 }
643
644 // last row elements, lacking a below row
645 for (int j = 0; j < nbCols - 1; ++j) {
646 preprocessed[k] = selector.select(elevations[k], elevations[k + 1]);
647 k++;
648 }
649
650 // last element
651 preprocessed[k] = elevations[k];
652
653 }
654
655 /** Recursive application of a function.
656 * @param tree tree to fill-up with the recursive applications
657 * @param level current level
658 * @param levelRows number of rows at current level
659 * @param levelColumns number of columns at current level
660 * @param selector to apply
661 * @param base base array from which function arguments are drawn
662 * @param first index of the first element to consider in base array
663 */
664 private void applyRecursively(final double[] tree,
665 final int level, final int levelRows, final int levelColumns,
666 final Selector selector,
667 final double[] base, final int first) {
668
669 if (isColumnMerging(level + 1)) {
670
671 // merge columns pairs
672 int iTree = start[level];
673 int iBase = first;
674 final int nextColumns = (levelColumns + 1) / 2;
675 final boolean odd = (levelColumns & 0x1) != 0;
676 final int jEnd = odd ? nextColumns - 1 : nextColumns;
677 for (int i = 0; i < levelRows; ++i) {
678
679 // regular pairs
680 for (int j = 0; j < jEnd; ++j) {
681 tree[iTree++] = selector.select(base[iBase], base[iBase + 1]);
682 iBase += 2;
683 }
684
685 if (odd) {
686 // last column
687 tree[iTree++] = base[iBase++];
688 }
689
690
691 }
692
693 if (level > 0) {
694 applyRecursively(tree, level - 1, levelRows, nextColumns, selector, tree, start[level]);
695 }
696
697 } else {
698
699 // merge rows pairs
700 int iTree = start[level];
701 int iBase = first;
702 final int nextRows = (levelRows + 1) / 2;
703 final boolean odd = (levelRows & 0x1) != 0;
704 final int iEnd = odd ? nextRows - 1 : nextRows;
705
706 // regular pairs
707 for (int i = 0; i < iEnd; ++i) {
708
709 for (int j = 0; j < levelColumns; ++j) {
710 tree[iTree++] = selector.select(base[iBase], base[iBase + levelColumns]);
711 iBase++;
712 }
713 iBase += levelColumns;
714
715 }
716
717 if (odd) {
718 // last row
719 System.arraycopy(base, iBase, tree, iTree, levelColumns);
720 }
721
722 if (level > 0) {
723 applyRecursively(tree, level - 1, nextRows, levelColumns, selector, tree, start[level]);
724 }
725
726 }
727 }
728
729 }