source: Revenant/geoworld/src/gwmorphs/winderode.cpp

#include "winderode.h"

#include <omp.h>

#include "geoworld.h"

namespace geoworld
{
        void WindErosionGM::init(SynthesisSession & manager, const IGeoWorldAccessor & accessor, const unsigned short nMaxWidth, const unsigned short nMaxHeight)
        {
                _psfDensity = accessor.queryField(FT_Density);
        }

        void WindErosionGM::doMorph( LockedGWSection & section ) const
        {
                using namespace mars;
                using namespace geoworld;

                LockedGWSection::View
                        * pView = section.createView(section.bbox.left, section.bbox.top, section.bbox.right, section.bbox.bottom);

                ScalarField< RationalType > 
                        hmf(pView->width, pView->height, Wrap),
                        hmdf(pView->width, pView->height, Wrap);

                const RangeX< GeoHeightMap::Precision > rnghm = pView->getHeightMap() ->range();
                GaussianFn< RationalType > 
                        gauss(1.0, 1.0, rnghm.delta * _frDecayLevel, 0);
                signed int i, j, k, g;
                VectorTag< RationalType >::V2 vgr;

                ptr< NoiseGen::NoiseMap > pnm = _ngen.compute(pView->width, pView->height);
                const NoiseGen::NoiseMap & nm = *pnm;
                const RangeX< NoiseGen::NoiseMap::Precision > rngnm = nm.range();
                const signed int len = static_cast< signed int > (pView->getHeightMap() ->length);      // It's ok, HeightMap.length is never more than the positive maximum of a signed-integer

                // Copy integral-precision field to rational-precision field
                hmf << *pView->getHeightMap();

                for (unsigned int c = 0; c < _iterations; ++c)
                {
                        for (j = 0; j < hmf.height; ++j)
                                for (i = 0; i < hmf.width; ++i)
                                {
                                        const RationalType mh = hmf.mean5(i, j);
                                        // This ensures that pits and gullies are ignored
                                        if (mh - hmf(i, j) < 0)
                                        {
                                                const RationalType      // Get the noise map ratio at i,j as a ratio to the noise map range
                                                        nr = static_cast< RationalType > ((nm(i,j) - rngnm.minimum) / rngnm.delta);

                                                // Obtain gradient into vgr, sampling distance of 5
                                                grad5(hmf, i, j, vgr); 

                                                const RationalType
                                                        // Sample density half-way down at current spot
                                                        fp = static_cast< RationalType > (1) / 
                                                                static_cast< RationalType > (
                                                                        (*_psfDensity) (i, j, static_cast< GeoHeightMap::Precision > (hmf(i,j) * 0.5f))
                                                                ),      
                                                        pgr = MAG(vgr),

                                                        // Amount of material to remove
                                                        rdh = pgr * fp * nr;

                                                // If the displacement is zero or the volume of material to move is zero then skip
                                                if (rdh > 0.0)
                                                {
                                                        // Remove a smooth pile of material from the source site
                                                        for (k = i - 2; k <= i + 2; ++k)
                                                                for (g = j - 2; g <= j + 2; ++g)
                                                                        hmdf(k, g) -= gaussblur(k - i, g - j) * rdh;

                                                        const RationalType
                                                                rddh = rdh * gauss.f(pgr);      // TODO: Export gauss fn to pre-init tables state

                                                        // Scale-up the displacement by a random amount, makes a very pleasing result
                                                        vgr *= RANDf(mh) + static_cast< RationalType > (1);

                                                        const signed int 
                                                                x0 = i - RNDi(vgr.x), 
                                                                y0 = j - RNDi(vgr.y);

                                                        // Deposits "alluvium" into a smooth pile
                                                        for (k = x0 - 2; k <= x0 + 2; ++k)
                                                                for (g = y0 - 2; g <= y0 + 2; ++g)
                                                                        hmdf(k, g) += gaussblur(k - x0, g - y0) * rddh;
                                                }
                                        }
                                }

                        hmf += hmdf;
                        hmdf.clear();
                }

                if (_bSmoothPost)
                {
                        #pragma omp parallel for private(i, j)
                        for (j = 0; j < hmf.height; ++j)
                                for (i = 0; i < hmf.width; ++i)
                                        (*pView->getHeightMap())(i, j) = static_cast< GeoHeightMap::Precision > (
                                                hmf(i,j) + (
                                                        static_cast< RationalType > (hmf.meanVN(i, j) - hmf(i, j))
                                                ) 
                                                * (nm(i, j) - rngnm.minimum) / rngnm.delta
                                        );
                } else
                        *pView->getHeightMap() << hmf;

                section.releaseView(pView);
        }

        WindErosionGM::~WindErosionGM()
        {
                // HACK: Use releaseField
                delete _psfDensity;
        }

        WindErosionGM::WindErosionGM( const unsigned int nIterations /*= 8*/, const float frDecayLevel /*= 0.1f*/, const bool bSmoothPost /*= true */ ) : _iterations(nIterations), _frDecayLevel(frDecayLevel), _ngen(0.3f, 0.5f), _bSmoothPost(bSmoothPost), _pAccessor(NULL), _psfDensity(NULL)
        {
                using namespace mars;
                GaussianFn< RationalType > gauss(
                        static_cast< RationalType > (1), 
                        static_cast< RationalType > (1) / sqrt (
                                static_cast< RationalType > (2) * static_cast< RationalType > (PI)
                        ), 
                        1, 
                        0
                );      // TODO: Export to pre-init tables API

                for (signed int x = -2; x <= +2; ++x)
                        for (signed int y = -2; y <= +2; ++y)
                                gaussblur(x, y) = gauss.f(static_cast< RationalType > (x), static_cast< RationalType > (y));
        }

        GeoMorph * WindErosionGTFactory::createRandomInstance() const
        { return new WindErosionGM(_mmnIterations.next(), _frDecayLevel, _bSmoothPost); }

        void WindErosionGTFactory::configure( IConfigGTFactory * pFactoryConfig, const IConfigGTFactory::Settings & settings )
        {
                const geoworld::IConfigSection & section = *pFactoryConfig->getSection();

                *section["iterations"] >> _mmnIterations;
                *section["decay"] >> _frDecayLevel;
                *section["smooth"] >> _bSmoothPost;
        }
}