source: Revenant/geoworld/include/lodtools.h

#ifndef __GWLODTOOLS_H__
#define __GWLODTOOLS_H__

#include <mars_grids.h>
#include "geoworld.h"

namespace geoworld
{
        template< typename T >
        class LODMagnifier
        {
        private:
                typedef mars::MidpointDisplacementGrid< T > SeedingMPDG;
        using NoiseFn = typename SeedingMPDG::BrownianFn;
                        
                SeedingMPDG 
                        ** _ppMPDG0,
                        * _pmpdg;
                TileBlock< T > _block;
                const mars::ScalarField< T > & _hmLowLOD;
                unsigned short 
                        _nDeltaLOD, 
                        _nTileDim;
                NoiseFn _functor;
                unsigned long
                        _nWWidth, _nWHeight;

                unsigned int
                        _nSectionTilesAcross, _nSectionTilesDown;

        public:
                LODMagnifier (const mars::ScalarField< T > & hmLowLOD, const unsigned short nDeltaLOD, const unsigned short nTileDim, const float fCoarseness, const unsigned long nWWidth, const unsigned long nWHeight)
                        : 
                        _ppMPDG0 (new SeedingMPDG *[nWWidth / hmLowLOD.width]),
                        _pmpdg (SeedingMPDG::createInstance(hmLowLOD.width, hmLowLOD.height, mars::Normal)),
                        _hmLowLOD(hmLowLOD), _nDeltaLOD(nDeltaLOD), _functor(NoiseFn(std::min(hmLowLOD.width, hmLowLOD.height), fCoarseness)),
                        _nTileDim(nTileDim), _nWWidth(nWWidth), _nWHeight(nWHeight), 
                        _nSectionTilesAcross(hmLowLOD.width / nTileDim), _nSectionTilesDown(hmLowLOD.height / nTileDim),
                        _block (nTileDim)
                {
                        assert(((nTileDim - 1) & nTileDim) == 0); // Must be a power of 2
                        assert(_nSectionTilesAcross * nTileDim == hmLowLOD.width && _nSectionTilesDown * nTileDim == hmLowLOD.height);
                        assert(_nWWidth >> nDeltaLOD == hmLowLOD.width && _nWHeight >> nDeltaLOD == hmLowLOD.height);

                        const unsigned int nSectionsWide = nWWidth / hmLowLOD.width;
                        for (unsigned int i = 0; i < nSectionsWide; ++i)
                                _ppMPDG0[i] = SeedingMPDG::createInstance(hmLowLOD.width, hmLowLOD.height, mars::Normal);
                }
                ~LODMagnifier()
                {
                        const unsigned int nSectionsWide = _nWWidth / _hmLowLOD.width;
                        for (unsigned int i = 0; i < nSectionsWide; ++i)
                                delete _ppMPDG0[i];

                        delete [] _ppMPDG0;
                        delete _pmpdg;
                }

                void write (FileSource::MainDEM & tmap)
                {
                        using namespace mars;

                        tmap.seekp(0, 0, std::ios_base::beg);

                        for (unsigned int jj = 0; jj < _nWHeight; jj += _hmLowLOD.height)
                                for (unsigned int ii = 0, c = 0; ii < _nWWidth; ii += _hmLowLOD.width, ++c)
                                {
                                        const unsigned int nStitch = (jj > 0 ? StitchTop : 0) | (ii > 0 ? StitchLeft : 0);
                                        SeedingMPDG * pMPDG = _ppMPDG0[c];

                                        if (jj > 0)
                                                pMPDG->stitch(StitchTop);
                                        if (c > 0)
                                                pMPDG->stitch(StitchLeft, _ppMPDG0[c - 1]);

                                        pMPDG->seedFrom (_hmLowLOD, ii >> _nDeltaLOD, jj >> _nDeltaLOD, _nDeltaLOD);

                                        pMPDG->reset();
                                        pMPDG->step(pMPDG->itercount - _nDeltaLOD);
                                        pMPDG->iterations (_nDeltaLOD, _functor, (jj > 0 ? StitchTop : 0) | (ii > 0 ? StitchLeft : 0));

                                        tmap.seekp(ii / _nTileDim, jj / _nTileDim, std::ios::beg);

                                        for (unsigned int j = 0; j < _nSectionTilesDown; )
                                        {
                                                for (unsigned int i = 0; i < _nSectionTilesAcross; ++i)
                                                {
                                                        _block.copyFrom(*pMPDG, i, j);
                                                        _block >> tmap;
                                                }
                                                if (++j < _nSectionTilesDown)
                                                        tmap.seekp(-static_cast< signed int > (_nSectionTilesAcross), 1, std::ios_base::cur);
                                        }
                                }
                }
        };
}

#endif