mercator  0.4.0
A terrain generation library for the Worldforge system.
HeightMap.cpp
1 // This file may be redistributed and modified only under the terms of
2 // the GNU General Public License (See COPYING for details).
3 // Copyright (C) 2003 Alistair Riddoch, Damien McGinnes
4 
5 #ifdef HAVE_CONFIG_H
6 #include "config.h"
7 #endif
8 
9 #include "iround.h"
10 
11 #include "HeightMap.h"
12 #include "Terrain.h"
13 #include "TerrainMod.h"
14 #include "Surface.h"
15 #include "BasePoint.h"
16 #include "Area.h"
17 #include "Shader.h"
18 
19 #include <wfmath/MersenneTwister.h>
20 
21 #include <cmath>
22 #include <cassert>
23 #include <algorithm>
24 
25 namespace Mercator {
26 
32 class LinInterp {
33  private:
35  float m_size;
37  bool noCalc;
38  public:
40  float ep1, ep2;
42  inline float calc(float loc)
43  {
44  return ((noCalc) ? ep1 :
45  ((m_size-loc) * ep1 + loc * ep2));
46  }
52  LinInterp(float size,float l, float h) : m_size(size), noCalc(false),
53  ep1(l/size), ep2(h/size)
54  {
55  if (l==h) {
56  ep1 = l;
57  noCalc=true;
58  }
59  }
60 };
61 
68 class QuadInterp {
69  private:
71  float m_size;
73  bool noCalc;
74  public:
76  float ep1, ep2, ep3, ep4;
78  inline float calc(float locX, float locY)
79  {
80  return ((noCalc) ? ep1 :
81  (( ep1*(m_size-locX) + ep2 * locX) * (m_size-locY) +
82  ( ep4*(m_size-locX) + ep3 * locX) * (locY) ) / m_size );
83  }
91  QuadInterp(float size,float e1, float e2, float e3, float e4)
92  : m_size(size), noCalc(false),
93  ep1(e1/size), ep2(e2/size), ep3(e3/size), ep4(e4/size)
94  {
95  if ((e1==e2) && (e3==e4) && (e2==e3)) {
96  ep1 = e1;
97  noCalc=true;
98  }
99  }
100 };
101 
102 
104 HeightMap::HeightMap(int resolution) :
105  Buffer<float>::Buffer(resolution + 1, 1),
106  m_res(resolution),
107  m_max(std::numeric_limits<float>::lowest()),
108  m_min(std::numeric_limits<float>::max())
109 {
110 }
111 
112 
118 {
119  if (h<m_min) {
120  m_min=h;
121  }
122  if (h>m_max) {
123  m_max=h;
124  }
125 }
126 
127 // generate a rand num between -0.5...0.5
128 inline float randHalf(WFMath::MTRand& rng)
129 {
130  //return (float) rand() / RAND_MAX - 0.5f;
131  return rng.rand<float>() - 0.5f;
132 }
133 
134 
136 float HeightMap::qRMD(WFMath::MTRand& rng, float nn, float fn, float ff, float nf,
137  float roughness, float falloff, float depth) const
138 {
139  float max = std::max(std::max(nn, fn), std::max(nf, ff)),
140  min = std::min(std::min(nn, fn), std::min(nf, ff)),
141  heightDifference = max - min;
142 
143  return ((nn+fn+ff+nf)/4.f) + randHalf(rng) * roughness * heightDifference / (1.f+std::pow(depth,falloff));
144 }
145 
152 void HeightMap::fill1d(const BasePoint& l, const BasePoint &h,
153  float *array) const
154 {
155  array[0] = l.height();
156  array[m_res] = h.height();
157  LinInterp li((float)m_res, l.roughness(), h.roughness());
158 
159  // seed the RNG.
160  // The RNG is seeded only once for the line and the seed is based on the
161  // two endpoints -because they are the common parameters for two adjoining
162  // tiles
163  //srand((l.seed() * 1000 + h.seed()));
164  WFMath::MTRand::uint32 seed[2]={ l.seed(), h.seed() };
165  WFMath::MTRand rng(seed, 2);
166 
167  // stride is used to step across the array in a deterministic fashion
168  // effectively we do the 1/2 point, then the 1/4 points, then the 1/8th
169  // points etc. this has to be the same order every time because we call
170  // on the RNG at every point
171  int stride = m_res/2;
172 
173  // depth is used to indicate what level we are on. the displacement is
174  // reduced each time we traverse the array.
175  float depth=1;
176 
177  while (stride) {
178  for (int i=stride;i<m_res;i+=stride*2) {
179  float hh = array[i-stride];
180  float lh = array[i+stride];
181  float hd = std::fabs(hh-lh);
182  float roughness = li.calc((float)i);
183 
184  //eliminate the problem where hd is nearly zero, leaving a flat section.
185  if ((hd*100.f) < roughness) {
186  hd+=0.05f * roughness;
187  }
188 
189  array[i] = ((hh+lh)/2.f) + randHalf(rng) * roughness * hd / (1.f+std::pow(depth,BasePoint::FALLOFF));
190  }
191  stride >>= 1;
192  depth++;
193  }
194 }
195 
201 void HeightMap::fill2d(const BasePoint& p1, const BasePoint& p2,
202  const BasePoint& p3, const BasePoint& p4)
203 {
204  //First reset the min and max values, since they will be updated.
205  m_max = std::numeric_limits<float>::lowest();
206  m_min = std::numeric_limits<float>::max();
207 
208  // calculate the edges first. This is necessary so that segments tile
209  // seamlessly note the order in which the edges are calculated and the
210  // direction. opposite edges are calculated the same way (eg left->right)
211  // so that the top of one tile matches the bottom of another, likewise
212  // with sides.
213 
214  // temporary array used to hold each edge
215  std::vector<float> edgeData;
216  edgeData.reserve(m_size);
217  float* edge = edgeData.data();
218 
219  float* points = m_data.data();
220 
221  // calc top edge and copy into m_heightMap
222  fill1d(p1,p2,edge);
223  for (int i=0;i<=m_res;i++) {
224  points[0*m_size + i] = edge[i];
225  checkMaxMin(edge[i]);
226  }
227 
228  // calc left edge and copy into points
229  fill1d(p1,p4,edge);
230  for (int i=0;i<=m_res;i++) {
231  points[i*m_size + 0] = edge[i];
232  checkMaxMin(edge[i]);
233  }
234 
235  // calc right edge and copy into points
236  fill1d(p2,p3,edge);
237  for (int i=0;i<=m_res;i++) {
238  points[i*m_size + m_res] = edge[i];
239  checkMaxMin(edge[i]);
240  }
241 
242  // calc bottom edge and copy into points
243  fill1d(p4,p3,edge);
244  for (int i=0;i<=m_res;i++) {
245  points[m_res*m_size + i] = edge[i];
246  checkMaxMin(edge[i]);
247  }
248 
249  // seed the RNG - this is the 5th and last seeding for the tile.
250  // it was seeded once for each edge, now once for the tile.
251  //srand(p1.seed()*20 + p2.seed()*15 + p3.seed()*10 + p4.seed()*5);
252  WFMath::MTRand::uint32 seed[4]={ p1.seed(), p2.seed(), p3.seed(), p4.seed() };
253  WFMath::MTRand rng(seed, 4);
254 
255  QuadInterp qi((float)m_res, p1.roughness(), p2.roughness(), p3.roughness(), p4.roughness());
256  QuadInterp falloffQi((float)m_res, p1.falloff(), p2.falloff(), p3.falloff(), p4.falloff());
257 
258  float depth=0;
259 
260  // center of points is done separately
261  int stride = m_res/2;
262 
263  //float roughness = (p1.roughness+p2.roughness+p3.roughness+p4.roughness)/(4.0f);
264  float roughness = qi.calc((float)stride,(float) stride);
265  float f = falloffQi.calc((float)stride, (float)stride);
266  points[stride*m_size + stride] = qRMD(rng, points[0 * m_size + stride],
267  points[stride*m_size + 0],
268  points[stride*m_size + m_res],
269  points[m_res*m_size + stride],
270  roughness,
271  f, depth);
272 
273 
274  checkMaxMin(points[stride*m_size + stride]);
275 
276  stride >>= 1;
277 
278  // skip across the points and fill in the points
279  // alternate cross and plus shapes.
280  // this is a diamond-square algorithm.
281  while (stride) {
282  //Cross shape - + contributes to value at X
283  //+ . +
284  //. X .
285  //+ . +
286  for (int i=stride;i<m_res;i+=stride*2) {
287  for (int j=stride;j<m_res;j+=stride*2) {
288  roughness=qi.calc((float)i,(float)j);
289  f = falloffQi.calc((float)i, (float)j);
290  points[j*m_size + i] = qRMD(rng, points[(i-stride) + (j+stride) * (m_size)],
291  points[(i+stride) + (j-stride) * (m_size)],
292  points[(i+stride) + (j+stride) * (m_size)],
293  points[(i-stride) + (j-stride) * (m_size)],
294  roughness, f, depth);
295  checkMaxMin(points[j*m_size + i]);
296  }
297  }
298 
299  depth++;
300  //Plus shape - + contributes to value at X
301  //. + .
302  //+ X +
303  //. + .
304  for (int i=stride*2;i<m_res;i+=stride*2) {
305  for (int j=stride;j<m_res;j+=stride*2) {
306  roughness=qi.calc((float)i,(float)j);
307  f = falloffQi.calc((float)i, (float)j);
308  points[j*m_size + i] = qRMD(rng, points[(i-stride) + (j) * (m_size)],
309  points[(i+stride) + (j) * (m_size)],
310  points[(i) + (j+stride) * (m_size)],
311  points[(i) + (j-stride) * (m_size)],
312  roughness, f , depth);
313  checkMaxMin(points[j*m_size + i]);
314  }
315  }
316 
317  for (int i=stride;i<m_res;i+=stride*2) {
318  for (int j=stride*2;j<m_res;j+=stride*2) {
319  roughness=qi.calc((float)i,(float)j);
320  f = falloffQi.calc((float)i, (float)j);
321  points[j*m_size + i] = qRMD(rng, points[(i-stride) + (j) * (m_size)],
322  points[(i+stride) + (j) * (m_size)],
323  points[(i) + (j+stride) * (m_size)],
324  points[(i) + (j-stride) * (m_size)],
325  roughness, f, depth);
326  checkMaxMin(points[j*m_size + i]);
327  }
328  }
329 
330  stride>>=1;
331  depth++;
332  }
333 }
334 
335 void HeightMap::getHeight(float x, float z, float &h) const
336 {
337  // FIXME this ignores edges and corners
338  assert(x <= m_res);
339  assert(x >= 0.0f);
340  assert(z <= m_res);
341  assert(z >= 0.0f);
342 
343  // get index of the actual tile in the segment
344  int tile_x = I_ROUND(std::floor(x));
345  int tile_z = I_ROUND(std::floor(z));
346 
347  // work out the offset into that tile
348  float off_x = x - (float)tile_x;
349  float off_z = z - (float)tile_z;
350 
351  float h1=get(tile_x, tile_z);
352  float h2=get(tile_x, tile_z+1);
353  float h3=get(tile_x+1, tile_z+1);
354  float h4=get(tile_x+1, tile_z);
355 
356  // square is broken into two triangles
357  // top triangle |/
358  if ((off_x - off_z) <= 0.f) {
359  h = h1 + (h3-h2) * off_x + (h2-h1) * off_z;
360  }
361  // bottom triangle /|
362  else {
363  h = h1 + (h4-h1) * off_x + (h3-h4) * off_z;
364  }
365 }
366 
367 
380 void HeightMap::getHeightAndNormal(float x, float z, float& h,
381  WFMath::Vector<3> &normal) const
382 {
383  // FIXME this ignores edges and corners
384  assert(x <= m_res);
385  assert(x >= 0.0f);
386  assert(z <= m_res);
387  assert(z >= 0.0f);
388 
389  // get index of the actual tile in the segment
390  int tile_x = I_ROUND(std::floor(x));
391  int tile_z = I_ROUND(std::floor(z));
392 
393  // work out the offset into that tile
394  float off_x = x - (float)tile_x;
395  float off_z = z - (float)tile_z;
396 
397  float h1=get(tile_x, tile_z);
398  float h2=get(tile_x, tile_z+1);
399  float h3=get(tile_x+1, tile_z+1);
400  float h4=get(tile_x+1, tile_z);
401 
402  // square is broken into two triangles
403  // top triangle |/
404  if ((off_x - off_z) <= 0.f) {
405  normal = WFMath::Vector<3>(h2-h3, 1.0f, h1-h2);
406 
407  //normal for intersection of both triangles
408  if (off_x == off_z) {
409  normal += WFMath::Vector<3>(h1-h4, 1.0f, h4-h3);
410  }
411  normal.normalize();
412  h = h1 + (h3-h2) * off_x + (h2-h1) * off_z;
413  }
414  // bottom triangle /|
415  else {
416  normal = WFMath::Vector<3>(h1-h4, 1.0f, h4-h3);
417  normal.normalize();
418  h = h1 + (h4-h1) * off_x + (h3-h4) * off_z;
419  }
420 }
421 
422 } // namespace Mercator
void fill2d(const BasePoint &p1, const BasePoint &p2, const BasePoint &p3, const BasePoint &p4)
Two dimensional midpoint displacement fractal.
Definition: HeightMap.cpp:201
void getHeightAndNormal(float x, float z, float &h, WFMath::Vector< 3 > &normal) const
Get an accurate height and normal vector at a given coordinate relative to this segment.
Definition: HeightMap.cpp:380
float ep1
Values at the two ends.
Definition: HeightMap.cpp:40
float calc(float loc)
Determine the interpolated value along the line.
Definition: HeightMap.cpp:42
Template for managing buffers of data for a segment.
Definition: Buffer.h:14
STL namespace.
QuadInterp(float size, float e1, float e2, float e3, float e4)
Constructor.
Definition: HeightMap.cpp:91
const unsigned int m_size
The size of segment, m_res + 1.
Definition: Buffer.h:19
Helper to interpolate in a quad.
Definition: HeightMap.cpp:68
LinInterp(float size, float l, float h)
Constructor.
Definition: HeightMap.cpp:52
float calc(float locX, float locY)
Determine the interpolated value within the quad.
Definition: HeightMap.cpp:78
Point on the fundamental grid that is used as the basis for terrain.
Definition: BasePoint.h:19
float roughness() const
Accessor for the roughness at the base point.
Definition: BasePoint.h:56
float ep1
Values at the four corners.
Definition: HeightMap.cpp:76
Helper to interpolate on a line.
Definition: HeightMap.cpp:32
std::vector< float > m_data
Pointer to buffer containing data values.
Definition: Buffer.h:21
static constexpr float FALLOFF
Default falloff at the base point.
Definition: BasePoint.h:34
float height() const
Accessor for the height at the base point.
Definition: BasePoint.h:51
void checkMaxMin(float h)
Check a value against m_min and m_max and set one of them if appropriate.
Definition: HeightMap.cpp:117
float falloff() const
Accessor for the falloff at the base point.
Definition: BasePoint.h:61
HeightMap(int resolution)
Construct an empty height map with the given resolution.
Definition: HeightMap.cpp:104
unsigned int seed() const
Calculate the random seed used at this base point.
Definition: BasePoint.cpp:14