I wrote an implementation, but the results are weird, can you help me find the problem? Below is a sample weird image (the sharp peaks and holes are troublesome to me) and code:
http://www.fasterupload.com/viewer.php?file=m8uumafcszrne1ox3hld.png
http://www.fasterupload.com/viewer.php?file=bqald61waxvrr3kbrxox.png
#ifndef MIDPOINTa_HPP
# define MIDPOINTa_HPP
#include <cmath>
#include "boost/multi_array.hpp"
#include "boost/random.hpp"
//////////////////////////////////////////////////////////////////////////////
template <class R, class T>
void diamond_stepa(R& vng, unsigned x, unsigned y, unsigned step,
boost::multi_array<T, 2>& hf)
{
T step2(step / 2);
T mdx(x + step2);
T mdy(y + step2);
hf[mdy][mdx] = (hf[y][x] + hf[y][x + step] + hf[y + step][x] +
hf[y + step][x + step]) / 4 + vng();
}
//////////////////////////////////////////////////////////////////////////////
template <class R, class T>
void square_stepa(R& vng, unsigned x, unsigned y, unsigned step,
boost::multi_array<T, 2>& hf)
{
T step2(step / 2);
T mdx(x + step2);
T mdy(y + step2);
// left
if (x)
{
hf[mdy][x] = (hf[y][x] + hf[y + step][x] + hf[mdy][mdx] +
hf[mdy][x - step2]) / 4 + vng();
}
else
{
hf[mdy][x] = (hf[y][x] + hf[y + step][x] + hf[mdy][mdx]) / 3 + vng();
}
// right
if (x + step == hf.size() - 1)
{
hf[mdy][x + step] = (hf[y][x + step] + hf[y + step][x + step] +
hf[mdy][mdx]) / 3 + vng();
}
// up
if (y)
{
hf[y][mdx] = (hf[y][x] + hf[y][x + step] + hf[mdy][mdx] +
hf[y - step2][mdx]) / 4 + vng();
}
else
{
hf[y][mdx] = (hf[y][x] + hf[y][x + step] + hf[mdy][mdx]) / 3 +
vng();
}
// down
if (y + step == hf.size() - 1)
{
hf[y + step][mdx] = (hf[y + step][x] + hf[y + step][x + step] +
hf[mdy][mdx]) / 3 + vng();
}
}
//////////////////////////////////////////////////////////////////////////////
template <class R, class T>
void midpointa(R& rng, unsigned n, T hval, T rc, boost::multi_array<T, 2>& hf)
{
hval /= 2;
boost::uniform_real<T> distrib(-hval, hval);
boost::variate_generator<R, boost::uniform_real<T> > vrg(rng, distrib);
unsigned step(unsigned(std::pow(2, n)) + 1);
typedef boost::multi_array<T , 2> array_type;
typename array_type::extent_gen extents;
hf.resize(extents[step][step]);
--step;
hf[0][0] = vrg();
hf[0][step] = vrg();
hf[step][0] = vrg();
hf[step][step] = vrg();
unsigned m(step);
do
{
for (unsigned x(0); x != m; x += step)
{
for (unsigned y(0); y != m; y += step)
{
diamond_stepa(vrg, x, y, step, hf);
}
}
for (unsigned x(0); x != m; x += step)
{
for (unsigned y(0); y != m; y += step)
{
square_stepa(vrg, x, y, step, hf);
}
}
hval *= std::pow(2, rc);
boost::uniform_real<T> distrib(-hval, hval);
vrg.distribution() = distrib;
step /= 2;
}
while (--n);
}
#endif // MIDPOINTa_HPP