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
25namespace Mercator {
26
32class 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
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
104HeightMap::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
128inline float randHalf(WFMath::MTRand& rng)
129{
130 //return (float) rand() / RAND_MAX - 0.5f;
131 return rng.rand<float>() - 0.5f;
132}
133
134
136float 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
152void 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
201void 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
335void 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
380void 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
Point on the fundamental grid that is used as the basis for terrain.
Definition: BasePoint.h:19
static constexpr float FALLOFF
Default falloff at the base point.
Definition: BasePoint.h:34
float falloff() const
Accessor for the falloff at the base point.
Definition: BasePoint.h:61
float roughness() const
Accessor for the roughness at the base point.
Definition: BasePoint.h:56
unsigned int seed() const
Calculate the random seed used at this base point.
Definition: BasePoint.cpp:14
Template for managing buffers of data for a segment.
Definition: Buffer.h:14
const unsigned int m_size
The size of segment, m_res + 1.
Definition: Buffer.h:19
std::vector< float > m_data
Pointer to buffer containing data values.
Definition: Buffer.h:21
void fill2d(const BasePoint &p1, const BasePoint &p2, const BasePoint &p3, const BasePoint &p4)
Two dimensional midpoint displacement fractal.
Definition: HeightMap.cpp:201
void checkMaxMin(float h)
Check a value against m_min and m_max and set one of them if appropriate.
Definition: HeightMap.cpp:117
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 get(int x, int z) const
Get the height at a relative integer position in the Segment.
Definition: HeightMap.h:40
HeightMap(int resolution)
Construct an empty height map with the given resolution.
Definition: HeightMap.cpp:104
Helper to interpolate on a line.
Definition: HeightMap.cpp:32
LinInterp(float size, float l, float h)
Constructor.
Definition: HeightMap.cpp:52
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
Helper to interpolate in a quad.
Definition: HeightMap.cpp:68
QuadInterp(float size, float e1, float e2, float e3, float e4)
Constructor.
Definition: HeightMap.cpp:91
float calc(float locX, float locY)
Determine the interpolated value within the quad.
Definition: HeightMap.cpp:78
float ep1
Values at the four corners.
Definition: HeightMap.cpp:76