text: /******************************************************************************
   *
   * 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 "cpl_port.h"
  #include "gdal_alg.h"
  
  #include <cmath>
  #include <cstring>
  #include <array>
  #include <limits>
  #include <algorithm>
  
  #include "cpl_conv.h"
  #include "cpl_error.h"
  #include "cpl_progress.h"
  #include "cpl_vsi.h"
  #include "gdal.h"
  #include "gdal_priv.h"
  #include "gdal_priv_templates.hpp"
  #include "ogr_api.h"
  #include "ogr_spatialref.h"
  #include "ogr_core.h"
  #include "commonutils.h"
  
  CPL_CVSID("$Id: viewshed.cpp 85dfaa122ba5861e5350f34850c217a8e52ae3f8 2021-12-26 10:08:17 +0100 Even Rouault $")
  
  
  inline static void SetVisibility(int iPixel, double dfZ, double dfZTarget, double* padfZVal,
      std::vector<GByte>& vResult, GByte byVisibleVal, GByte byInvisibleVal)
  {
      if (padfZVal[iPixel] + dfZTarget < dfZ)
          vResult[iPixel] = byInvisibleVal;
      else
          vResult[iPixel] = byVisibleVal;
  
      if (padfZVal[iPixel] < dfZ)
          padfZVal[iPixel] = dfZ;
  }
  
  inline static bool AdjustHeightInRange(const double* adfGeoTransform, int iPixel, int iLine, double& dfHeight, double dfDistance2, double dfCurvCoeff, double dfSphereDiameter)
  {
      if (dfDistance2 <= 0 && dfCurvCoeff == 0)
          return true;
  
      double dfX = adfGeoTransform[1] * iPixel + adfGeoTransform[2] * iLine;
      double dfY = adfGeoTransform[4] * iPixel + adfGeoTransform[5] * iLine;
      double dfR2 = dfX * dfX + dfY * dfY;
  
      /* calc adjustment */
      if (dfCurvCoeff != 0 && dfSphereDiameter != std::numeric_limits<double>::infinity())
          dfHeight -= dfCurvCoeff * dfR2 / dfSphereDiameter;
  
      if (dfDistance2 > 0 && dfR2 > dfDistance2)
          return false;
  
      return true;
  }
  
  inline static double CalcHeightLine(int i, double Za, double Zo)
  {
      if (i == 1)
          return Za;
      else
          return (Za - Zo) / (i - 1) + Za;
  }
  
  inline static double CalcHeightDiagonal(int i, int j, double Za, double Zb, double Zo)
  {
      return ((Za - Zo) * i + (Zb - Zo) * j) / (i + j - 1) + Zo;
  }
  
  inline static double CalcHeightEdge(int i, int j, double Za, double Zb, double Zo)
  {
      if (i == j)
          return CalcHeightLine(i, Za, Zo);
      else
          return ((Za - Zo) * i + (Zb - Zo) * (j - i)) / (j - 1) + Zo;
  }
  
  inline static double CalcHeight(double dfZ, double dfZ2, GDALViewshedMode eMode)
  {
      if (eMode == GVM_Edge)
          return dfZ2;
      else if (eMode == GVM_Max)
          return (std::max)(dfZ, dfZ2);
      else if (eMode == GVM_Min)
          return (std::min)(dfZ, dfZ2);
      else
          return dfZ;
  }
  
  
  /************************************************************************/
  /*                        GDALViewshedGenerate()                         */
  /************************************************************************/
  
  /**
   * Create viewshed from raster DEM.
   *
   * This algorithm will generate a viewshed raster from an input DEM raster
   * by using a modified algorithm of "Generating Viewsheds without Using Sightlines"
   * published at https://www.asprs.org/wp-content/uploads/pers/2000journal/january/2000_jan_87-90.pdf
   * This appoach provides a relatively fast calculation, since the output raster is
   * generated in a single scan.
   * The gdal/apps/gdal_viewshed.cpp mainline can be used as an example of
   * how to use this function.
   * The output raster will be of type Byte or Float64.
   *
   * \note The algorithm as implemented currently will only output meaningful results
   * if the georeferencing is in a projected coordinate reference system.
   *
   * @param hBand The band to read the DEM data from. Only the part of the raster
   * within the specified maxdistance around the observer point is processed.
   *
   * @param pszDriverName Driver name (GTiff if set to NULL)
   *
   * @param pszTargetRasterName The name of the target raster to be generated. Must not be NULL
   *
   * @param papszCreationOptions creation options.
   *
   * @param dfObserverX observer X value (in SRS units)
   *
   * @param dfObserverY observer Y value (in SRS units)
   *
   * @param dfObserverHeight The height of the observer above the DEM surface.
   *
   * @param dfTargetHeight The height of the target above the DEM surface. (default 0)
   *
   * @param dfVisibleVal pixel value for visibility (default 255)
   *
   * @param dfInvisibleVal pixel value for invisibility (default 0)
   *
   * @param dfOutOfRangeVal The value to be set for the cells that fall outside of the
   * range specified by dfMaxDistance.
   *
   * @param dfNoDataVal The value to be set for the cells that have no data.
   *                    If set to a negative value, nodata is not set.
   *                    Note: currently, no special processing of input cells at a nodata
   *                    value is done (which may result in erroneous results).
   *
   * @param dfCurvCoeff Coefficient to consider the effect of the curvature and refraction.
   * The height of the DEM is corrected according to the following formula:
   * [Height] -= dfCurvCoeff * [Target Distance]^2 / [Earth Diameter]
   * For the effect of the atmospheric refraction we can use 0.85714.
   *
   * @param eMode The mode of the viewshed calculation.
   * Possible values GVM_Diagonal = 1, GVM_Edge = 2 (default), GVM_Max = 3, GVM_Min = 4.
   *
   * @param dfMaxDistance maximum distance range to compute viewshed.
   *                      It is also used to clamp the extent of the output raster.
   *                      If set to 0, then unlimited range is assumed, that is to say the
   *                      computation is performed on the extent of the whole raster.
   *
   * @param pfnProgress A GDALProgressFunc that may be used to report progress
   * to the user, or to interrupt the algorithm.  May be NULL if not required.
   *
   * @param pProgressArg The callback data for the pfnProgress function.
   *
   * @param heightMode Type of information contained in output raster. Possible values
   *                   GVOT_NORMAL = 1 (default), GVOT_MIN_TARGET_HEIGHT_FROM_DEM = 2,
   *                   GVOT_MIN_TARGET_HEIGHT_FROM_GROUND = 3
   *
   *                   GVOT_NORMAL returns a raster of type Byte containing visible locations.
   *
   *                   GVOT_MIN_TARGET_HEIGHT_FROM_DEM and GVOT_MIN_TARGET_HEIGHT_FROM_GROUND
   *                   will return a raster of type Float64 containing the minimum target height
   *                   for target to be visible from the DEM surface or ground level respectively.
   *                   Parameters dfTargetHeight, dfVisibleVal and dfInvisibleVal will be ignored.
   *
   *
   * @param papszExtraOptions Future extra options. Must be set to NULL currently.
   *
   * @return not NULL output dataset on success (to be closed with GDALClose()) or NULL if an error occurs.
   *
   * @since GDAL 3.1
   */
  
  GDALDatasetH GDALViewshedGenerate(GDALRasterBandH hBand,
                              const char* pszDriverName,
                              const char* pszTargetRasterName,
                              CSLConstList papszCreationOptions,
                      double dfObserverX, double dfObserverY, double dfObserverHeight,
                      double dfTargetHeight, double dfVisibleVal, double dfInvisibleVal,
                      double dfOutOfRangeVal, double dfNoDataVal, double dfCurvCoeff,
                      GDALViewshedMode eMode, double dfMaxDistance,
                      GDALProgressFunc pfnProgress, void *pProgressArg,
                      GDALViewshedOutputType heightMode, CSLConstList papszExtraOptions)
  
  {
      VALIDATE_POINTER1( hBand, "GDALViewshedGenerate", nullptr );
      VALIDATE_POINTER1( pszTargetRasterName, "GDALViewshedGenerate", nullptr );
  
      CPL_IGNORE_RET_VAL(papszExtraOptions);
  
      if( pfnProgress == nullptr )
          pfnProgress = GDALDummyProgress;
  
      if( !pfnProgress( 0.0, "", pProgressArg ) )
      {
          CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
          return nullptr;
      }
  
      const GByte byNoDataVal = dfNoDataVal >= 0 && dfNoDataVal <= 255 ? static_cast<GByte>(dfNoDataVal) : 0;
      const GByte byVisibleVal = dfVisibleVal >= 0 && dfVisibleVal <= 255 ? static_cast<GByte>(dfVisibleVal) : 255;
      const GByte byInvisibleVal = dfInvisibleVal >= 0 && dfInvisibleVal <= 255 ? static_cast<GByte>(dfInvisibleVal) : 0;
      const GByte byOutOfRangeVal = dfOutOfRangeVal >= 0 && dfOutOfRangeVal <= 255 ? static_cast<GByte>(dfOutOfRangeVal) : 0;
  
      if(heightMode != GVOT_MIN_TARGET_HEIGHT_FROM_DEM && heightMode != GVOT_MIN_TARGET_HEIGHT_FROM_GROUND)
          heightMode = GVOT_NORMAL;
  
      /* set up geotransformation */
      std::array<double, 6> adfGeoTransform {{0.0, 1.0, 0.0, 0.0, 0.0, 1.0}};
      GDALDatasetH hSrcDS = GDALGetBandDataset( hBand );
      if( hSrcDS != nullptr )
          GDALGetGeoTransform( hSrcDS, adfGeoTransform.data());
  
      double adfInvGeoTransform[6];
      if (!GDALInvGeoTransform(adfGeoTransform.data(), adfInvGeoTransform))
      {
          CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
          return nullptr;
      }
  
      /* calculate observer position */
      double dfX, dfY;
      GDALApplyGeoTransform(adfInvGeoTransform, dfObserverX, dfObserverY, &dfX, &dfY);
      int nX = static_cast<int>(dfX);
      int nY = static_cast<int>(dfY);
  
      int nXSize = GDALGetRasterBandXSize( hBand );
      int nYSize = GDALGetRasterBandYSize( hBand );
  
      if (nX < 0 ||
          nX > nXSize ||
          nY < 0 ||
          nY > nYSize)
      {
          CPLError(CE_Failure, CPLE_AppDefined, "The observer location falls outside of the DEM area");
          return nullptr;
      }
  
      /* calculate the area of interest */
      int nXStart = dfMaxDistance > 0? (std::max)(0, static_cast<int>(std::floor(nX - adfInvGeoTransform[1] * dfMaxDistance))) : 0;
      int nXStop = dfMaxDistance > 0? (std::min)(nXSize, static_cast<int>(std::ceil(nX + adfInvGeoTransform[1] * dfMaxDistance) + 1)) : nXSize;
      int nYStart = dfMaxDistance > 0? (std::max)(0, static_cast<int>(std::floor(nY + adfInvGeoTransform[5] * dfMaxDistance))) : 0;
      int nYStop = dfMaxDistance > 0? (std::min)(nYSize, static_cast<int>(std::ceil(nY - adfInvGeoTransform[5] * dfMaxDistance) + 1)) : nYSize;
  
      /* normalize horizontal index (0 - nXSize) */
      nXSize = nXStop - nXStart;
      nX -= nXStart;
  
      nYSize = nYStop - nYStart;
  
      if (nXSize == 0 || nYSize == 0)
      {
          CPLError(CE_Failure, CPLE_AppDefined, "Invalid target raster size");
          return nullptr;
      }
  
      std::vector<double> vFirstLineVal;
      std::vector<double> vLastLineVal;
      std::vector<double> vThisLineVal;
      std::vector<GByte> vResult;
      std::vector<double> vHeightResult;
  
      try
      {
          vFirstLineVal.resize(nXSize);
          vLastLineVal.resize(nXSize);
          vThisLineVal.resize(nXSize);
          vResult.resize(nXSize);
  
          if(heightMode != GVOT_NORMAL)
              vHeightResult.resize(nXSize);
      } catch (...)
      {
          CPLError(CE_Failure, CPLE_AppDefined,
                   "Cannot allocate vectors for viewshed");
          return nullptr;
      }
  
      double *padfFirstLineVal = vFirstLineVal.data();
      double *padfLastLineVal = vLastLineVal.data();
      double *padfThisLineVal = vThisLineVal.data();
      GByte *pabyResult = vResult.data();
      double *dfHeightResult = vHeightResult.data();
  
      GDALDriverManager *hMgr = GetGDALDriverManager();
      GDALDriver *hDriver = hMgr->GetDriverByName(pszDriverName ? pszDriverName : "GTiff");
      if (!hDriver)
      {
          CPLError(CE_Failure, CPLE_AppDefined, "Cannot get driver");
          return nullptr;
      }
  
      /* create output raster */
      auto poDstDS = std::unique_ptr<GDALDataset>(hDriver->Create(pszTargetRasterName, nXSize, nYStop - nYStart, 1, heightMode != GVOT_NORMAL ? GDT_Float64 : GDT_Byte,
                                                  const_cast<char**>(papszCreationOptions)));
      if (!poDstDS)
      {
          CPLError(CE_Failure, CPLE_AppDefined,
              "Cannot create dataset for %s", pszTargetRasterName);
          return nullptr;
      }
      /* copy srs */
      if (hSrcDS)
          poDstDS->SetSpatialRef(GDALDataset::FromHandle(hSrcDS)->GetSpatialRef());
  
      std::array<double, 6> adfDstGeoTransform;
      adfDstGeoTransform[0] = adfGeoTransform[0] + adfGeoTransform[1] * nXStart + adfGeoTransform[2] * nYStart;
      adfDstGeoTransform[1] = adfGeoTransform[1];
      adfDstGeoTransform[2] = adfGeoTransform[2];
      adfDstGeoTransform[3] = adfGeoTransform[3] + adfGeoTransform[4] * nXStart + adfGeoTransform[5] * nYStart;
      adfDstGeoTransform[4] = adfGeoTransform[4];
      adfDstGeoTransform[5] = adfGeoTransform[5];
      poDstDS->SetGeoTransform(adfDstGeoTransform.data());
  
      auto hTargetBand = poDstDS->GetRasterBand(1);
      if (hTargetBand == nullptr)
      {
          CPLError(CE_Failure, CPLE_AppDefined,
              "Cannot get band for %s", pszTargetRasterName);
          return nullptr;
      }
  
      if (dfNoDataVal >= 0)
          GDALSetRasterNoDataValue(hTargetBand, heightMode != GVOT_NORMAL ? dfNoDataVal : byNoDataVal);
  
      /* process first line */
      if (GDALRasterIO(hBand, GF_Read, nXStart, nY, nXSize, 1,
          padfFirstLineVal, nXSize, 1, GDT_Float64, 0, 0))
      {
          CPLError(CE_Failure, CPLE_AppDefined,
              "RasterIO error when reading DEM at position(%d, %d), size(%d, %d)", nXStart, nY, nXSize, 1);
          return nullptr;
      }
  
      const double dfZObserver = dfObserverHeight + padfFirstLineVal[nX];
      double dfZ = 0.0;
      const double dfDistance2 = dfMaxDistance * dfMaxDistance;
  
      /* If we can't get a SemiMajor axis from the SRS, it will be
       * SRS_WGS84_SEMIMAJOR
      */
      double dfSphereDiameter(std::numeric_limits<double>::infinity());
      const OGRSpatialReference* poDstSRS = poDstDS->GetSpatialRef();
      if (poDstSRS)
      {
          OGRErr eSRSerr;
          double dfSemiMajor = poDstSRS->GetSemiMajor(&eSRSerr);
  
          /* If we fetched the axis from the SRS, use it */
          if (eSRSerr != OGRERR_FAILURE)
              dfSphereDiameter = dfSemiMajor * 2.0;
          else
              CPLDebug( "GDALViewshedGenerate", "Unable to fetch SemiMajor axis from spatial reference");
  
      }
  
      /* mark the observer point as visible */
      double dfGroundLevel = heightMode == GVOT_MIN_TARGET_HEIGHT_FROM_DEM ? padfFirstLineVal[nX] : 0.0;
      pabyResult[nX] = byVisibleVal;
      if(heightMode != GVOT_NORMAL)
          dfHeightResult[nX] = dfGroundLevel;
  
      if (nX > 0)
      {
          dfGroundLevel = heightMode == GVOT_MIN_TARGET_HEIGHT_FROM_DEM ? padfFirstLineVal[nX - 1] : 0.0;
          CPL_IGNORE_RET_VAL(
              AdjustHeightInRange(adfGeoTransform.data(),
                              1,
                              0,
                              padfFirstLineVal[nX - 1],
                              dfDistance2,
                              dfCurvCoeff,
                              dfSphereDiameter));
          pabyResult[nX - 1] = byVisibleVal;
          if(heightMode != GVOT_NORMAL)
              dfHeightResult[nX - 1] = dfGroundLevel;
      }
      if (nX < nXSize - 1)
      {
          dfGroundLevel = heightMode == GVOT_MIN_TARGET_HEIGHT_FROM_DEM ? padfFirstLineVal[nX + 1] : 0.0;
          CPL_IGNORE_RET_VAL(
              AdjustHeightInRange(adfGeoTransform.data(),
                              1,
                              0,
                              padfFirstLineVal[nX + 1],
                              dfDistance2,
                              dfCurvCoeff,
                              dfSphereDiameter));
          pabyResult[nX + 1] = byVisibleVal;
          if(heightMode != GVOT_NORMAL)
              dfHeightResult[nX + 1] = dfGroundLevel;
      }
  
      /* process left direction */
      for (int iPixel = nX - 2; iPixel >= 0; iPixel--)
      {
          dfGroundLevel = heightMode == GVOT_MIN_TARGET_HEIGHT_FROM_DEM ? padfFirstLineVal[iPixel] : 0.0;
          bool adjusted = AdjustHeightInRange(adfGeoTransform.data(),
                                              nX - iPixel,
                                              0,
                                              padfFirstLineVal[iPixel],
                                              dfDistance2,
                                              dfCurvCoeff,
                                              dfSphereDiameter);
          if (adjusted)
          {
              dfZ = CalcHeightLine(nX - iPixel,
                                   padfFirstLineVal[iPixel + 1],
                                   dfZObserver);
  
              if(heightMode != GVOT_NORMAL)
                  dfHeightResult[iPixel] = std::max(0.0, (dfZ - padfFirstLineVal[iPixel] + dfGroundLevel));
  
              SetVisibility(  iPixel,
                              dfZ,
                              dfTargetHeight,
                              padfFirstLineVal,
                              vResult,
                              byVisibleVal,
                              byInvisibleVal);
          }
          else
          {
              for (; iPixel >= 0; iPixel--)
              {
                  pabyResult[iPixel] = byOutOfRangeVal;
                  if(heightMode != GVOT_NORMAL)
                      dfHeightResult[iPixel] = dfOutOfRangeVal;
              }
          }
      }
      /* process right direction */
      for (int iPixel = nX + 2; iPixel < nXSize; iPixel++)
      {
          dfGroundLevel = heightMode == GVOT_MIN_TARGET_HEIGHT_FROM_DEM ? padfFirstLineVal[iPixel] : 0.0;
          bool adjusted = AdjustHeightInRange(adfGeoTransform.data(),
                                              iPixel - nX,
                                              0,
                                              padfFirstLineVal[iPixel],
                                              dfDistance2,
                                              dfCurvCoeff,
                                              dfSphereDiameter);
          if (adjusted)
          {
              dfZ = CalcHeightLine(iPixel - nX,
                                   padfFirstLineVal[iPixel - 1],
                                   dfZObserver);
  
              if(heightMode != GVOT_NORMAL)
                  dfHeightResult[iPixel] = std::max(0.0, (dfZ - padfFirstLineVal[iPixel] + dfGroundLevel));
  
              SetVisibility(iPixel,
                            dfZ,
                            dfTargetHeight,
                            padfFirstLineVal,
                            vResult,
                            byVisibleVal,
                            byInvisibleVal);
          }
          else
          {
              for (; iPixel < nXSize; iPixel++)
              {
                  pabyResult[iPixel] = byOutOfRangeVal;
                  if(heightMode != GVOT_NORMAL)
                      dfHeightResult[iPixel] = dfOutOfRangeVal;
              }
          }
      }
      /* write result line */
  
      if (GDALRasterIO(hTargetBand, GF_Write, 0, nY - nYStart, nXSize, 1,
          heightMode != GVOT_NORMAL ? static_cast<void*>(dfHeightResult) : static_cast<void*>(pabyResult), nXSize, 1, heightMode != GVOT_NORMAL ? GDT_Float64 : GDT_Byte, 0, 0))
      {
          CPLError(CE_Failure, CPLE_AppDefined,
              "RasterIO error when writing target raster at position (%d,%d), size (%d,%d)", 0, nY - nYStart, nXSize, 1);
          return nullptr;
      }
  
      /* scan upwards */
      std::copy(vFirstLineVal.begin(),
                vFirstLineVal.end(),
                vLastLineVal.begin());
      for (int iLine = nY - 1; iLine >= nYStart; iLine--)
      {
          if (GDALRasterIO(hBand, GF_Read, nXStart, iLine, nXSize, 1,
              padfThisLineVal, nXSize, 1, GDT_Float64, 0, 0))
          {
              CPLError(CE_Failure, CPLE_AppDefined,
                  "RasterIO error when reading DEM at position (%d,%d), size (%d,%d)", nXStart, iLine, nXSize, 1);
              return nullptr;
          }
  
          /* set up initial point on the scanline */
          dfGroundLevel = heightMode == GVOT_MIN_TARGET_HEIGHT_FROM_DEM ? padfThisLineVal[nX] : 0.0;
          bool adjusted = AdjustHeightInRange(adfGeoTransform.data(),
                                              0,
                                              nY - iLine,
                                              padfThisLineVal[nX],
                                              dfDistance2,
                                              dfCurvCoeff,
                                              dfSphereDiameter);
          if (adjusted)
          {
              dfZ = CalcHeightLine(nY - iLine,
                                   padfLastLineVal[nX],
                                   dfZObserver);
  
              if(heightMode != GVOT_NORMAL)
                  dfHeightResult[nX] = std::max(0.0, (dfZ - padfThisLineVal[nX] + dfGroundLevel));
  
              SetVisibility(nX,
                            dfZ,
                            dfTargetHeight,
                            padfThisLineVal,
                            vResult,
                            byVisibleVal,
                            byInvisibleVal);
          }
          else
          {
              pabyResult[nX] = byOutOfRangeVal;
              if(heightMode != GVOT_NORMAL)
                  dfHeightResult[nX] = dfOutOfRangeVal;
          }
  
          /* process left direction */
          for (int iPixel = nX - 1; iPixel >= 0; iPixel--)
          {
              dfGroundLevel = heightMode == GVOT_MIN_TARGET_HEIGHT_FROM_DEM ? padfThisLineVal[iPixel] : 0.0;
              bool left_adjusted = AdjustHeightInRange(adfGeoTransform.data(),
                                                       nX - iPixel,
                                                       nY - iLine,
                                                       padfThisLineVal[iPixel],
                                                       dfDistance2,
                                                       dfCurvCoeff,
                                                       dfSphereDiameter);
              if (left_adjusted)
              {
                  if (eMode != GVM_Edge)
                      dfZ = CalcHeightDiagonal(nX - iPixel,
                                               nY - iLine,
                                               padfThisLineVal[iPixel + 1],
                                               padfLastLineVal[iPixel],
                                               dfZObserver);
  
                  if (eMode != GVM_Diagonal)
                  {
                      double dfZ2 = nX - iPixel >= nY - iLine ?
                          CalcHeightEdge(nY - iLine,
                                         nX - iPixel,
                                         padfLastLineVal[iPixel + 1],
                                         padfThisLineVal[iPixel + 1],
                                         dfZObserver) :
                          CalcHeightEdge(nX - iPixel,
                                         nY - iLine,
                                         padfLastLineVal[iPixel + 1],
                                         padfLastLineVal[iPixel],
                                         dfZObserver);
                      dfZ = CalcHeight(dfZ, dfZ2, eMode);
                  }
  
                  if(heightMode != GVOT_NORMAL)
                      dfHeightResult[iPixel] = std::max(0.0, (dfZ - padfThisLineVal[iPixel] + dfGroundLevel));
  
                  SetVisibility(iPixel,
                                dfZ,
                                dfTargetHeight,
                                padfThisLineVal,
                                vResult,
                                byVisibleVal,
                                byInvisibleVal);
              }
              else
              {
                  for (; iPixel >= 0; iPixel--)
                  {
                      pabyResult[iPixel] = byOutOfRangeVal;
                      if(heightMode != GVOT_NORMAL)
                          dfHeightResult[iPixel] = dfOutOfRangeVal;
                  }
              }
          }
          /* process right direction */
          for (int iPixel = nX + 1; iPixel < nXSize; iPixel++)
          {
              dfGroundLevel = heightMode == GVOT_MIN_TARGET_HEIGHT_FROM_DEM ? padfThisLineVal[iPixel] : 0.0;
              bool right_adjusted = AdjustHeightInRange(adfGeoTransform.data(),
                                                        iPixel - nX,
                                                        nY - iLine,
                                                        padfThisLineVal[iPixel],
                                                        dfDistance2,
                                                        dfCurvCoeff,
                                                        dfSphereDiameter);
              if (right_adjusted)
              {
                  if (eMode != GVM_Edge)
                      dfZ = CalcHeightDiagonal(iPixel - nX,
                                               nY - iLine,
                                               padfThisLineVal[iPixel - 1],
                                               padfLastLineVal[iPixel],
                                               dfZObserver);
  
                  if (eMode != GVM_Diagonal)
                  {
                      double dfZ2 = iPixel - nX >= nY - iLine ?
                          CalcHeightEdge(nY - iLine,
                                         iPixel - nX,
                                         padfLastLineVal[iPixel - 1],
                                         padfThisLineVal[iPixel - 1],
                                         dfZObserver) :
                          CalcHeightEdge(iPixel - nX,
                                         nY - iLine,
                                         padfLastLineVal[iPixel - 1],
                                         padfLastLineVal[iPixel],
                                         dfZObserver);
                      dfZ = CalcHeight(dfZ, dfZ2, eMode);
                  }
  
                  if(heightMode != GVOT_NORMAL)
                      dfHeightResult[iPixel] = std::max(0.0, (dfZ - padfThisLineVal[iPixel] + dfGroundLevel));
  
                  SetVisibility(iPixel,
                                dfZ,
                                dfTargetHeight,
                                padfThisLineVal,
                                vResult,
                                byVisibleVal,
                                byInvisibleVal);
              }
              else
              {
                  for (; iPixel < nXSize; iPixel++)
                  {
                      pabyResult[iPixel] = byOutOfRangeVal;
                      if(heightMode != GVOT_NORMAL)
                          dfHeightResult[iPixel] = dfOutOfRangeVal;
                  }
              }
          }
  
          /* write result line */
          if (GDALRasterIO(hTargetBand, GF_Write, 0, iLine - nYStart, nXSize, 1,
              heightMode != GVOT_NORMAL ? static_cast<void*>(dfHeightResult) : static_cast<void*>(pabyResult), nXSize, 1, heightMode != GVOT_NORMAL ? GDT_Float64 : GDT_Byte, 0, 0))
          {
              CPLError(CE_Failure, CPLE_AppDefined,
                  "RasterIO error when writing target raster at position (%d,%d), size (%d,%d)", 0, iLine - nYStart, nXSize, 1);
              return nullptr;
          }
  
          std::swap(padfLastLineVal, padfThisLineVal);
  
          if (!pfnProgress((nY - iLine) / static_cast<double>(nYSize),
                  "", pProgressArg))
          {
              CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
              return nullptr;
          }
      }
      /* scan downwards */
      memcpy(padfLastLineVal, padfFirstLineVal, nXSize * sizeof(double));
      for(int iLine = nY + 1; iLine < nYStop; iLine++ )
      {
          if (GDALRasterIO( hBand, GF_Read, nXStart, iLine, nXStop - nXStart, 1,
              padfThisLineVal, nXStop - nXStart, 1, GDT_Float64, 0, 0 ))
          {
              CPLError(CE_Failure, CPLE_AppDefined,
                  "RasterIO error when reading DEM at position (%d,%d), size (%d,%d)", nXStart, iLine, nXStop - nXStart, 1);
              return nullptr;
          }
  
          /* set up initial point on the scanline */
          dfGroundLevel = heightMode == GVOT_MIN_TARGET_HEIGHT_FROM_DEM ? padfThisLineVal[nX] : 0.0;
          bool adjusted = AdjustHeightInRange(adfGeoTransform.data(),
                                              0,
                                              iLine - nY,
                                              padfThisLineVal[nX],
                                              dfDistance2,
                                              dfCurvCoeff,
                                              dfSphereDiameter);
          if (adjusted)
          {
              dfZ = CalcHeightLine(iLine - nY,
                                   padfLastLineVal[nX],
                                   dfZObserver);
  
              if(heightMode != GVOT_NORMAL)
                  dfHeightResult[nX] = std::max(0.0, (dfZ - padfThisLineVal[nX] + dfGroundLevel));
  
              SetVisibility(nX,
                            dfZ,
                            dfTargetHeight,
                            padfThisLineVal,
                            vResult,
                            byVisibleVal,
                            byInvisibleVal);
          }
          else
          {
              pabyResult[nX] = byOutOfRangeVal;
              if(heightMode != GVOT_NORMAL)
                  dfHeightResult[nX] = dfOutOfRangeVal;
          }
  
          /* process left direction */
          for (int iPixel = nX - 1; iPixel >= 0; iPixel--)
          {
              dfGroundLevel = heightMode == GVOT_MIN_TARGET_HEIGHT_FROM_DEM ? padfThisLineVal[iPixel] : 0.0;
              bool left_adjusted = AdjustHeightInRange(adfGeoTransform.data(),
                                                       nX - iPixel,
                                                       iLine - nY,
                                                       padfThisLineVal[iPixel],
                                                       dfDistance2,
                                                       dfCurvCoeff,
                                                       dfSphereDiameter);
              if (left_adjusted)
              {
                  if (eMode != GVM_Edge)
                      dfZ = CalcHeightDiagonal(nX - iPixel, iLine - nY,
                          padfThisLineVal[iPixel + 1], padfLastLineVal[iPixel], dfZObserver);
  
                  if (eMode != GVM_Diagonal)
                  {
                      double dfZ2 = nX - iPixel >= iLine - nY ?
                          CalcHeightEdge(iLine - nY,
                                         nX - iPixel,
                                         padfLastLineVal[iPixel + 1],
                                         padfThisLineVal[iPixel + 1],
                                         dfZObserver) :
                          CalcHeightEdge(nX - iPixel,
                                         iLine - nY,
                                         padfLastLineVal[iPixel + 1],
                                         padfLastLineVal[iPixel],
                                         dfZObserver);
                      dfZ = CalcHeight(dfZ, dfZ2, eMode);
                  }
  
                  if(heightMode != GVOT_NORMAL)
                      dfHeightResult[iPixel] = std::max(0.0, (dfZ - padfThisLineVal[iPixel] + dfGroundLevel));
  
                  SetVisibility(iPixel,
                                dfZ,
                                dfTargetHeight,
                                padfThisLineVal,
                                vResult,
                                byVisibleVal,
                                byInvisibleVal);
              }
              else
              {
                  for (; iPixel >= 0; iPixel--)
                  {
                      pabyResult[iPixel] = byOutOfRangeVal;
                      if(heightMode != GVOT_NORMAL)
                          dfHeightResult[iPixel] = dfOutOfRangeVal;
                  }
              }
          }
          /* process right direction */
          for (int iPixel = nX + 1; iPixel < nXSize; iPixel++)
          {
              dfGroundLevel = heightMode == GVOT_MIN_TARGET_HEIGHT_FROM_DEM ? padfThisLineVal[iPixel] : 0.0;
              bool right_adjusted = AdjustHeightInRange(adfGeoTransform.data(),
                                                        iPixel - nX,
                                                        iLine - nY,
                                                        padfThisLineVal[iPixel],
                                                        dfDistance2,
                                                        dfCurvCoeff,
                                                        dfSphereDiameter);
              if (right_adjusted)
              {
                  if (eMode != GVM_Edge)
                      dfZ = CalcHeightDiagonal(iPixel - nX,
                                               iLine - nY,
                                               padfThisLineVal[iPixel - 1],
                                               padfLastLineVal[iPixel],
                                               dfZObserver);
  
                  if (eMode != GVM_Diagonal)
                  {
                      double dfZ2 = iPixel - nX >= iLine - nY ?
                          CalcHeightEdge(iLine - nY,
                                         iPixel - nX,
                                         padfLastLineVal[iPixel - 1],
                                         padfThisLineVal[iPixel - 1],
                                         dfZObserver) :
                          CalcHeightEdge(iPixel - nX,
                                         iLine - nY,
                                         padfLastLineVal[iPixel - 1],
                                         padfLastLineVal[iPixel],
                                         dfZObserver);
                      dfZ = CalcHeight(dfZ, dfZ2, eMode);
                  }
  
                  if(heightMode != GVOT_NORMAL)
                      dfHeightResult[iPixel] = std::max(0.0, (dfZ - padfThisLineVal[iPixel] + dfGroundLevel));
  
                  SetVisibility(iPixel,
                                dfZ,
                                dfTargetHeight,
                                padfThisLineVal,
                                vResult,
                                byVisibleVal,
                                byInvisibleVal);
              }
              else
              {
                  for (; iPixel < nXSize; iPixel++)
                  {
                      pabyResult[iPixel] = byOutOfRangeVal;
                      if(heightMode != GVOT_NORMAL)
                          dfHeightResult[iPixel] = dfOutOfRangeVal;
                  }
              }
          }
  
          /* write result line */
          if (GDALRasterIO(hTargetBand, GF_Write, 0, iLine - nYStart, nXSize, 1,
              heightMode != GVOT_NORMAL ? static_cast<void*>(dfHeightResult) : static_cast<void*>(pabyResult), nXSize, 1, heightMode != GVOT_NORMAL ? GDT_Float64 : GDT_Byte, 0, 0))
          {
              CPLError(CE_Failure, CPLE_AppDefined,
                  "RasterIO error when writing target raster at position (%d,%d), size (%d,%d)", 0, iLine - nYStart, nXSize, 1);
              return nullptr;
          }
  
          std::swap(padfLastLineVal, padfThisLineVal);
  
          if (!pfnProgress((iLine - nYStart) / static_cast<double>(nYSize),
              "", pProgressArg))
          {
              CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
              return nullptr;
          }
      }
  
      if (!pfnProgress(1.0, "", pProgressArg))
      {
          CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
          return nullptr;
      }
  
      return GDALDataset::FromHandle(poDstDS.release());
  }
  