viewshed.h

/******************************************************************************
 *
 * Project:  Viewshed Generation
 * Purpose:  Core algorithm implementation for viewshed generation.
 * Author:   Tamas Szekeres, szekerest@gmail.com
 *
 ******************************************************************************
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 ****************************************************************************/
 
#include <algorithm>
#include <array>
#include <cstdint>
#include <functional>
#include <limits>
#include <memory>
#include <mutex>
#include <string>
 
#include "cpl_progress.h"
#include "gdal_priv.h"
 
namespace gdal
{
 
class Viewshed
{
  public:
    enum class OutputMode
    {
        Normal,  
        DEM,     
        Ground   
    };
 
    enum class CellMode
    {
        Diagonal,  
        Edge,      
        Max,       
        Min        
    };
 
    struct Point
    {
        double x;  
        double y;  
        double z;  
    };
 
    struct Window
    {
        int xStart{};  
        int xStop{};   
        int yStart{};  
        int yStop{};   
 
        int xSize() const
        {
            return xStop - xStart;
        }
 
        int ySize() const
        {
            return yStop - yStart;
        }
 
        bool containsX(int nX) const
        {
            return nX >= xStart && nX < xStop;
        }
 
        bool containsY(int nY) const
        {
            return nY >= xStart && nY < yStop;
        }
 
        bool contains(int nX, int nY) const
        {
            return containsX(nX) && containsY(nY);
        }
 
        int clampX(int nX) const
        {
            return xSize() ? std::clamp(nX, xStart, xStop - 1) : xStart;
        }
 
        int clampY(int nY) const
        {
            return ySize() ? std::clamp(nY, yStart, yStop - 1) : yStart;
        }
 
        void shiftX(int nShift)
        {
            xStart += nShift;
            xStop += nShift;
        }
    };
 
    struct Options
    {
        Point observer{0, 0, 0};  
        double visibleVal{255};   
        double invisibleVal{
            0};  
        double outOfRangeVal{
            0};  
        double nodataVal{-1};  
        double targetHeight{0.0};  
        double maxDistance{
            0.0};  
        double curveCoeff{.85714};  
        OutputMode outputMode{OutputMode::Normal};  
        std::string outputFormat{};    
        std::string outputFilename{};  
        CPLStringList creationOpts{};  
        CellMode cellMode{
            CellMode::Edge};  
    };
 
    CPL_DLL explicit Viewshed(const Options &opts)
        : oOpts{opts}, oOutExtent{}, oCurExtent{},
          dfMaxDistance2{opts.maxDistance * opts.maxDistance},
          dfZObserver{0}, poDstDS{}, pSrcBand{}, pDstBand{},
          dfHeightAdjFactor{0}, nLineCount{0}, adfTransform{0, 1, 0, 0, 0, 1},
          adfInvTransform{}, oProgress{}, oZcalc{}, oMutex{}, iMutex{}
    {
        if (dfMaxDistance2 == 0)
            dfMaxDistance2 = std::numeric_limits<double>::max();
    }
 
    Viewshed(const Viewshed &) = delete;
    Viewshed &operator=(const Viewshed &) = delete;
 
    CPL_DLL bool run(GDALRasterBandH hBand,
                     GDALProgressFunc pfnProgress = GDALDummyProgress,
                     void *pProgressArg = nullptr);
 
    CPL_DLL std::unique_ptr<GDALDataset> output()
    {
        return std::move(poDstDS);
    }
 
  private:
    Options oOpts;
    Window oOutExtent;
    Window oCurExtent;
    double dfMaxDistance2;
    double dfZObserver;
    std::unique_ptr<GDALDataset> poDstDS;
    GDALRasterBand *pSrcBand;
    GDALRasterBand *pDstBand;
    double dfHeightAdjFactor;
    int nLineCount;
    std::array<double, 6> adfTransform;
    std::array<double, 6> adfInvTransform;
    using ProgressFunc = std::function<bool(double frac, const char *msg)>;
    ProgressFunc oProgress;
    using ZCalc = std::function<double(int, int, double, double, double)>;
    ZCalc oZcalc;
    std::mutex oMutex;
    std::mutex iMutex;
 
    void setOutput(double &dfResult, double &dfCellVal, double dfZ);
    double calcHeight(double dfZ, double dfZ2);
    bool readLine(int nLine, double *data);
    bool writeLine(int nLine, std::vector<double> &vResult);
    bool processLine(int nX, int nY, int nLine,
                     std::vector<double> &vLastLineVal);
    bool processFirstLine(int nX, int nY, std::vector<double> &vLastLineVal);
    void processFirstLineLeft(int nX, int iStart, int iEnd,
                              std::vector<double> &vResult,
                              std::vector<double> &vThisLineVal);
    void processFirstLineRight(int nX, int iStart, int iEnd,
                               std::vector<double> &vResult,
                               std::vector<double> &vThisLineVal);
    void processFirstLineTopOrBottom(int iLeft, int iRight,
                                     std::vector<double> &vResult,
                                     std::vector<double> &vThisLineVal);
    void processLineLeft(int nX, int nYOffset, int iStart, int iEnd,
                         std::vector<double> &vResult,
                         std::vector<double> &vThisLineVal,
                         std::vector<double> &vLastLineVal);
    void processLineRight(int nX, int nYOffset, int iStart, int iEnd,
                          std::vector<double> &vResult,
                          std::vector<double> &vThisLineVal,
                          std::vector<double> &vLastLineVal);
    std::pair<int, int> adjustHeight(int iLine, int nX,
                                     std::vector<double> &thisLineVal);
    bool calcOutputExtent(int nX, int nY);
    bool createOutputDataset();
    bool lineProgress();
    bool emitProgress(double fraction);
};
 
}  // namespace gdal