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.raster;
18  
19  import org.hipparchus.exception.LocalizedCoreFormats;
20  import org.hipparchus.exception.MathIllegalArgumentException;
21  import org.hipparchus.random.RandomGenerator;
22  import org.hipparchus.random.Well19937a;
23  import org.hipparchus.util.ArithmeticUtils;
24  import org.hipparchus.util.FastMath;
25  
26  public class RandomLandscapeUpdater implements TileUpdater {
27  
28      private double size;
29      private int    n;
30      private double[][] heightMap;
31  
32      /**
33       * @param baseH elevation of the base of DEM; unit = m
34       * @param initialScale
35       * @param reductionFactor for smoothing for low value (0.1) or not (0.5 for instance) the landscape
36       * @param seed
37       * @param size size in latitude / size in longitude (rad)
38       * @param n number of latitude / number of longitude
39       */
40      public RandomLandscapeUpdater(double baseH, double initialScale, double reductionFactor,
41                                    long seed, double size, int n) {
42  
43          if (!ArithmeticUtils.isPowerOfTwo(n - 1)) {
44              throw new MathIllegalArgumentException(LocalizedCoreFormats.SIMPLE_MESSAGE,
45                              "tile size must be a power of two plus one");
46          }
47  
48          this.size = size;
49          this.n    = n;
50  
51          // as we want to use this for testing and comparison purposes,
52          // and as we don't control when tiles are generated, we generate
53          // only ONCE a height map with continuity at the edges, and
54          // reuse this single height map for ALL generated tiles
55          heightMap = new double[n][n];
56          RandomGenerator random  = new Well19937a(seed);
57  
58          // initialize corners in diamond-square algorithm
59          heightMap[0][0]         = baseH;
60          heightMap[0][n - 1]     = baseH;
61          heightMap[n - 1][0]     = baseH;
62          heightMap[n - 1][n - 1] = baseH;
63  
64          double scale = initialScale;
65          for (int span = n - 1; span > 1; span = span / 2) {
66  
67              // perform squares step
68              for (int i = span / 2; i < n; i += span) {
69                  for (int j = span / 2; j < n; j += span) {
70                      double middleH = mean(i - span / 2, j - span / 2,
71                                            i - span / 2, j + span / 2,
72                                            i + span / 2, j - span / 2,
73                                            i + span / 2, j + span / 2) +
74                                      scale * (random.nextDouble() - 0.5);
75                      heightMap[i][j] = middleH;
76                  }
77              }
78  
79              // perform diamonds step
80              boolean flipFlop = false;
81              for (int i = 0; i < n - 1; i += span / 2) {
82                  for (int j = flipFlop ? 0 : span / 2; j < n - 1; j += span) {
83                      double middleH = mean(i - span / 2, j,
84                                            i + span / 2, j,
85                                            i,            j - span / 2,
86                                            i,            j + span / 2) +
87                                      scale * (random.nextDouble() - 0.5);
88                      heightMap[i][j] = middleH;
89                      if (i == 0) {
90                          heightMap[n - 1][j] = middleH;
91                      }
92                      if (j == 0) {
93                          heightMap[i][n - 1] = middleH;
94                      }
95                  }
96                  flipFlop = !flipFlop;
97              }
98  
99              // reduce scale
100             scale *= reductionFactor;
101 
102         }
103 
104     }
105 
106     @Override
107     public void updateTile(double latitude, double longitude, UpdatableTile tile) {
108                 
109         double step         = size / (n - 1);
110         double minLatitude  = size * FastMath.floor(latitude  / size);
111         double minLongitude = size * FastMath.floor(longitude / size);
112         tile.setGeometry(minLatitude, minLongitude, step, step, n, n);
113         for (int i = 0; i < n; ++i) {
114             for (int j = 0; j < n; ++j) {
115                 tile.setElevation(i, j, heightMap[i][j]);
116             }
117         }
118     }
119 
120     private double mean(int i1, int j1, int i2, int j2, int i3, int j3, int i4, int j4) {
121         return (heightMap[(i1 + n) % n][(j1 + n) % n] +
122                         heightMap[(i2 + n) % n][(j2 + n) % n] +
123                         heightMap[(i3 + n) % n][(j3 + n) % n] +
124                         heightMap[(i4 + n) % n][(j4 + n) % n]) / 4;
125     }
126 
127 }