15885686134

Committed 25 Jun 2025 07:44PM UTC coverage: 71.084%. Remained the same

Build # 15885686134

Build Type

push

github

Committed by

rouault

Commit Message

gdal_priv.h: fix C++11 compatibility

Run Details

573814 of 807237 relevant lines covered (71.08%)

250621.56 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

56.18

/alg/viewshed/viewshed.cpp

/******************************************************************************
 *
 * Project:  Viewshed Generation
 * Purpose:  Core algorithm implementation for viewshed generation.
 * Author:   Tamas Szekeres, szekerest@gmail.com
 *
 * (c) 2024 info@hobu.co
 *
 ******************************************************************************
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include <algorithm>
#include <array>

#include "gdal_alg.h"
#include "gdal_priv_templates.hpp"

#include "progress.h"
#include "util.h"
#include "viewshed.h"
#include "viewshed_executor.h"

/************************************************************************/
/*                        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 approach 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,
    [[maybe_unused]] CSLConstList papszExtraOptions)
{
    using namespace gdal;

    viewshed::Options oOpts;
    oOpts.outputFormat = pszDriverName;
    oOpts.outputFilename = pszTargetRasterName;
    oOpts.creationOpts = papszCreationOptions;
    oOpts.observer.x = dfObserverX;
    oOpts.observer.y = dfObserverY;
    oOpts.observer.z = dfObserverHeight;
    oOpts.targetHeight = dfTargetHeight;
    oOpts.curveCoeff = dfCurvCoeff;
    oOpts.maxDistance = dfMaxDistance;
    oOpts.nodataVal = dfNoDataVal;

    switch (eMode)
    {
        case GVM_Edge:
            oOpts.cellMode = viewshed::CellMode::Edge;
            break;
        case GVM_Diagonal:
            oOpts.cellMode = viewshed::CellMode::Diagonal;
            break;
        case GVM_Min:
            oOpts.cellMode = viewshed::CellMode::Min;
            break;
        case GVM_Max:
            oOpts.cellMode = viewshed::CellMode::Max;
            break;
    }

    switch (heightMode)
    {
        case GVOT_MIN_TARGET_HEIGHT_FROM_DEM:
            oOpts.outputMode = viewshed::OutputMode::DEM;
            break;
        case GVOT_MIN_TARGET_HEIGHT_FROM_GROUND:
            oOpts.outputMode = viewshed::OutputMode::Ground;
            break;
        case GVOT_NORMAL:
            oOpts.outputMode = viewshed::OutputMode::Normal;
            break;
    }

    if (!GDALIsValueInRange<uint8_t>(dfVisibleVal))
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "dfVisibleVal out of range. Must be [0, 255].");
        return nullptr;
    }
    if (!GDALIsValueInRange<uint8_t>(dfInvisibleVal))
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "dfInvisibleVal out of range. Must be [0, 255].");
        return nullptr;
    }
    if (!GDALIsValueInRange<uint8_t>(dfOutOfRangeVal))
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "dfOutOfRangeVal out of range. Must be [0, 255].");
        return nullptr;
    }
    oOpts.visibleVal = dfVisibleVal;
    oOpts.invisibleVal = dfInvisibleVal;
    oOpts.outOfRangeVal = dfOutOfRangeVal;

    gdal::viewshed::Viewshed v(oOpts);

    if (!pfnProgress)
        pfnProgress = GDALDummyProgress;
    v.run(hBand, pfnProgress, pProgressArg);

    return GDALDataset::FromHandle(v.output().release());
}

namespace gdal
{
namespace viewshed
{

namespace
{

bool getTransforms(GDALRasterBand &band, GDALGeoTransform &fwdTransform,
                   GDALGeoTransform &revTransform)
{
    band.GetDataset()->GetGeoTransform(fwdTransform);
    if (!fwdTransform.GetInverse(revTransform))
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
        return false;
    }
    return true;
}

/// Shrink the extent of a window to just cover the slice defined by rays from
/// (nX, nY) and [startAngle, endAngle]
///
/// @param oOutExtent  Window to modify
/// @param nX  X coordinate of ray endpoint.
/// @param nY  Y coordinate of ray endpoint.
/// @param startAngle  Start angle of slice (standard mathmatics notion, in radians)
/// @param endAngle  End angle of slice (standard mathmatics notion, in radians)
void shrinkWindowForAngles(Window &oOutExtent, int nX, int nY,
                           double startAngle, double endAngle)
{
    /// NOTE: This probably doesn't work when the observer is outside the raster and
    ///   needs to be enhanced for that case.

    if (startAngle == endAngle)
        return;

    Window win = oOutExtent;

    // Set the X boundaries for the angles
    int startAngleX = hIntersect(startAngle, nX, nY, win);
    int stopAngleX = hIntersect(endAngle, nX, nY, win);

    int xmax = nX;
    if (!rayBetween(startAngle, endAngle, 0))
    {
        xmax = std::max(xmax, startAngleX);
        xmax = std::max(xmax, stopAngleX);
        // Add one to xmax since we want one past the stop. [start, stop)
        oOutExtent.xStop = std::min(oOutExtent.xStop, xmax + 1);
    }

    int xmin = nX;
    if (!rayBetween(startAngle, endAngle, M_PI))
    {
        xmin = std::min(xmin, startAngleX);
        xmin = std::min(xmin, stopAngleX);
        oOutExtent.xStart = std::max(oOutExtent.xStart, xmin);
    }

    // Set the Y boundaries for the angles
    int startAngleY = vIntersect(startAngle, nX, nY, win);
    int stopAngleY = vIntersect(endAngle, nX, nY, win);

    int ymin = nY;
    if (!rayBetween(startAngle, endAngle, M_PI / 2))
    {
        ymin = std::min(ymin, startAngleY);
        ymin = std::min(ymin, stopAngleY);
        oOutExtent.yStart = std::max(oOutExtent.yStart, ymin);
    }
    int ymax = nY;
    if (!rayBetween(startAngle, endAngle, 3 * M_PI / 2))
    {
        ymax = std::max(ymax, startAngleY);
        ymax = std::max(ymax, stopAngleY);
        // Add one to ymax since we want one past the stop. [start, stop)
        oOutExtent.yStop = std::min(oOutExtent.yStop, ymax + 1);
    }
}

}  // unnamed namespace

Viewshed::Viewshed(const Options &opts) : oOpts{opts}
{
}

Viewshed::~Viewshed() = default;

/// Calculate the extent of the output raster in terms of the input raster and
/// save the input raster extent.
///
/// @return  false on error, true otherwise
bool Viewshed::calcExtents(int nX, int nY, const GDALGeoTransform &invGT)
{
    // We start with the assumption that the output size matches the input.
    oOutExtent.xStop = GDALGetRasterBandXSize(pSrcBand);
    oOutExtent.yStop = GDALGetRasterBandYSize(pSrcBand);

    if (!oOutExtent.contains(nX, nY))
    {
        if (oOpts.startAngle != oOpts.endAngle)
        {
            CPLError(CE_Failure, CPLE_AppDefined,
                     "Angle masking is not supported with an out-of-raster "
                     "observer.");
            return false;
        }
        CPLError(CE_Warning, CPLE_AppDefined,
                 "NOTE: The observer location falls outside of the DEM area");
    }

    constexpr double EPSILON = 1e-8;
    if (oOpts.maxDistance > 0)
    {
        //ABELL - This assumes that the transformation is only a scaling. Should be fixed.
        //  Find the distance in the direction of the transformed unit vector in the X and Y
        //  directions and use those factors to determine the limiting values in the raster space.
        int nXStart = static_cast<int>(
            std::floor(nX - invGT[1] * oOpts.maxDistance + EPSILON));
        int nXStop = static_cast<int>(
            std::ceil(nX + invGT[1] * oOpts.maxDistance - EPSILON) + 1);
        //ABELL - These seem to be wrong. The transform of 1 is no transform, so not
        //  sure why we're adding one in the first case. Really, the transformed distance
        // should add EPSILON. Not sure what the change should be for a negative transform,
        // which is what I think is being handled with the 1/0 addition/subtraction.
        int nYStart =
            static_cast<int>(std::floor(
                nY - std::fabs(invGT[5]) * oOpts.maxDistance + EPSILON)) -
            (invGT[5] > 0 ? 1 : 0);
        int nYStop = static_cast<int>(
            std::ceil(nY + std::fabs(invGT[5]) * oOpts.maxDistance - EPSILON) +
            (invGT[5] < 0 ? 1 : 0));

        // If the limits are invalid, set the window size to zero to trigger the error below.
        if (nXStart >= oOutExtent.xStop || nXStop < 0 ||
            nYStart >= oOutExtent.yStop || nYStop < 0)
        {
            oOutExtent = Window();
        }
        else
        {
            oOutExtent.xStart = std::max(nXStart, 0);
            oOutExtent.xStop = std::min(nXStop, oOutExtent.xStop);

            oOutExtent.yStart = std::max(nYStart, 0);
            oOutExtent.yStop = std::min(nYStop, oOutExtent.yStop);
        }
    }

    if (oOutExtent.xSize() == 0 || oOutExtent.ySize() == 0)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Invalid target raster size due to transform "
                 "and/or distance limitation.");
        return false;
    }

    shrinkWindowForAngles(oOutExtent, nX, nY, oOpts.startAngle, oOpts.endAngle);

    // normalize horizontal index to [ 0, oOutExtent.xSize() )
    oCurExtent = oOutExtent;
    oCurExtent.shiftX(-oOutExtent.xStart);

    return true;
}

/// Compute the viewshed of a raster band.
///
/// @param band  Pointer to the raster band to be processed.
/// @param pfnProgress  Pointer to the progress function. Can be null.
/// @param pProgressArg  Argument passed to the progress function
/// @return  True on success, false otherwise.
bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,
                   void *pProgressArg)
{
    pSrcBand = static_cast<GDALRasterBand *>(band);

    GDALGeoTransform fwdTransform, invTransform;
    if (!getTransforms(*pSrcBand, fwdTransform, invTransform))
        return false;

    // calculate observer position
    double dfX, dfY;
    invTransform.Apply(oOpts.observer.x, oOpts.observer.y, &dfX, &dfY);
    if (!GDALIsValueInRange<int>(dfX))
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Observer X value out of range");
        return false;
    }
    if (!GDALIsValueInRange<int>(dfY))
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Observer Y value out of range");
        return false;
    }
    int nX = static_cast<int>(dfX);
    int nY = static_cast<int>(dfY);

    if (oOpts.startAngle < 0 || oOpts.startAngle >= 360)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Start angle out of range. Must be [0, 360).");
        return false;
    }
    if (oOpts.endAngle < 0 || oOpts.endAngle >= 360)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "End angle out of range. Must be [0, 360).");
        return false;
    }
    if (oOpts.highPitch > 90)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Invalid highPitch. Cannot be greater than 90.");
        return false;
    }
    if (oOpts.lowPitch < -90)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Invalid lowPitch. Cannot be less than -90.");
        return false;
    }
    if (oOpts.highPitch <= oOpts.lowPitch)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Invalid pitch. highPitch must be > lowPitch");
        return false;
    }

    // Normalize angle to radians and standard math arrangement.
    oOpts.startAngle = normalizeAngle(oOpts.startAngle);
    oOpts.endAngle = normalizeAngle(oOpts.endAngle);

    // Must calculate extents in order to make the output dataset.
    if (!calcExtents(nX, nY, invTransform))
        return false;

    poDstDS = createOutputDataset(*pSrcBand, oOpts, oOutExtent);
    if (!poDstDS)
        return false;

    // Create the progress reporter.
    Progress oProgress(pfnProgress, pProgressArg, oOutExtent.ySize());

    // Execute the viewshed algorithm.
    GDALRasterBand *pDstBand = poDstDS->GetRasterBand(1);
    ViewshedExecutor executor(*pSrcBand, *pDstBand, nX, nY, oOutExtent,
                              oCurExtent, oOpts, oProgress,
                              /* emitWarningIfNoData = */ true);
    executor.run();
    oProgress.emit(1);
    return static_cast<bool>(poDstDS);
}

// Adjust the coefficient of curvature for non-earth SRS.
/// \param curveCoeff  Current curve coefficient
/// \param hSrcDS  Source dataset
/// \return  Adjusted curve coefficient.
double adjustCurveCoeff(double curveCoeff, GDALDatasetH hSrcDS)
{
    const OGRSpatialReference *poSRS =
        GDALDataset::FromHandle(hSrcDS)->GetSpatialRef();
    if (poSRS)
    {
        OGRErr eSRSerr;
        const double dfSemiMajor = poSRS->GetSemiMajor(&eSRSerr);
        if (eSRSerr != OGRERR_FAILURE &&
            fabs(dfSemiMajor - SRS_WGS84_SEMIMAJOR) >
                0.05 * SRS_WGS84_SEMIMAJOR)
        {
            curveCoeff = 1.0;
            CPLDebug("gdal_viewshed",
                     "Using -cc=1.0 as a non-Earth CRS has been detected");
        }
    }
    return curveCoeff;
}

#if defined(__clang__) || defined(__GNUC__)
#pragma GCC diagnostic ignored "-Wmissing-declarations"
#endif

void testShrinkWindowForAngles(Window &oOutExtent, int nX, int nY,
                               double startAngle, double endAngle)
{
    shrinkWindowForAngles(oOutExtent, nX, nY, startAngle, endAngle);
}

}  // namespace viewshed
}  // namespace gdal

1	/******************************************************************************
2	*
3	* Project: Viewshed Generation
4	* Purpose: Core algorithm implementation for viewshed generation.
5	* Author: Tamas Szekeres, szekerest@gmail.com
6	*
7	* (c) 2024 info@hobu.co
8	*
9	******************************************************************************
10	*
11	* SPDX-License-Identifier: MIT
12	****************************************************************************/
13
14	#include <algorithm>
15	#include <array>
16
17	#include "gdal_alg.h"
18	#include "gdal_priv_templates.hpp"
19
20	#include "progress.h"
21	#include "util.h"
22	#include "viewshed.h"
23	#include "viewshed_executor.h"
24
25	/************************************************************************/
26	/* GDALViewshedGenerate() */
27	/************************************************************************/
28
29	/**
30	* Create viewshed from raster DEM.
31	*
32	* This algorithm will generate a viewshed raster from an input DEM raster
33	* by using a modified algorithm of "Generating Viewsheds without Using
34	* Sightlines" published at
35	* https://www.asprs.org/wp-content/uploads/pers/2000journal/january/2000_jan_87-90.pdf
36	* This approach provides a relatively fast calculation, since the output raster
37	* is generated in a single scan. The gdal/apps/gdal_viewshed.cpp mainline can
38	* be used as an example of how to use this function. The output raster will be
39	* of type Byte or Float64.
40	*
41	* \note The algorithm as implemented currently will only output meaningful
42	* results if the georeferencing is in a projected coordinate reference system.
43	*
44	* @param hBand The band to read the DEM data from. Only the part of the raster
45	* within the specified maxdistance around the observer point is processed.
46	*
47	* @param pszDriverName Driver name (GTiff if set to NULL)
48	*
49	* @param pszTargetRasterName The name of the target raster to be generated.
50	* Must not be NULL
51	*
52	* @param papszCreationOptions creation options.
53	*
54	* @param dfObserverX observer X value (in SRS units)
55	*
56	* @param dfObserverY observer Y value (in SRS units)
57	*
58	* @param dfObserverHeight The height of the observer above the DEM surface.
59	*
60	* @param dfTargetHeight The height of the target above the DEM surface.
61	* (default 0)
62	*
63	* @param dfVisibleVal pixel value for visibility (default 255)
64	*
65	* @param dfInvisibleVal pixel value for invisibility (default 0)
66	*
67	* @param dfOutOfRangeVal The value to be set for the cells that fall outside of
68	* the range specified by dfMaxDistance.
69	*
70	* @param dfNoDataVal The value to be set for the cells that have no data.
71	* If set to a negative value, nodata is not set.
72	* Note: currently, no special processing of input cells at a
73	* nodata value is done (which may result in erroneous results).
74	*
75	* @param dfCurvCoeff Coefficient to consider the effect of the curvature and
76	* refraction. The height of the DEM is corrected according to the following
77	* formula: [Height] -= dfCurvCoeff * [Target Distance]^2 / [Earth Diameter] For
78	* the effect of the atmospheric refraction we can use 0.85714.
79	*
80	* @param eMode The mode of the viewshed calculation.
81	* Possible values GVM_Diagonal = 1, GVM_Edge = 2 (default), GVM_Max = 3,
82	* GVM_Min = 4.
83	*
84	* @param dfMaxDistance maximum distance range to compute viewshed.
85	* It is also used to clamp the extent of the output
86	* raster. If set to 0, then unlimited range is assumed, that is to say the
87	* computation is performed on the extent of the whole
88	* raster.
89	*
90	* @param pfnProgress A GDALProgressFunc that may be used to report progress
91	* to the user, or to interrupt the algorithm. May be NULL if not required.
92	*
93	* @param pProgressArg The callback data for the pfnProgress function.
94	*
95	* @param heightMode Type of information contained in output raster. Possible
96	* values GVOT_NORMAL = 1 (default), GVOT_MIN_TARGET_HEIGHT_FROM_DEM = 2,
97	* GVOT_MIN_TARGET_HEIGHT_FROM_GROUND = 3
98	*
99	* GVOT_NORMAL returns a raster of type Byte containing
100	* visible locations.
101	*
102	* GVOT_MIN_TARGET_HEIGHT_FROM_DEM and
103	* GVOT_MIN_TARGET_HEIGHT_FROM_GROUND will return a raster of type Float64
104	* containing the minimum target height for target to be visible from the DEM
105	* surface or ground level respectively. Parameters dfTargetHeight, dfVisibleVal
106	* and dfInvisibleVal will be ignored.
107	*
108	*
109	* @param papszExtraOptions Future extra options. Must be set to NULL currently.
110	*
111	* @return not NULL output dataset on success (to be closed with GDALClose()) or
112	* NULL if an error occurs.
113	*
114	* @since GDAL 3.1
115	*/
116	GDALDatasetH GDALViewshedGenerate(	×
117	GDALRasterBandH hBand, const char *pszDriverName,
118	const char *pszTargetRasterName, CSLConstList papszCreationOptions,
119	double dfObserverX, double dfObserverY, double dfObserverHeight,
120	double dfTargetHeight, double dfVisibleVal, double dfInvisibleVal,
121	double dfOutOfRangeVal, double dfNoDataVal, double dfCurvCoeff,
122	GDALViewshedMode eMode, double dfMaxDistance, GDALProgressFunc pfnProgress,
123	void *pProgressArg, GDALViewshedOutputType heightMode,
124	[[maybe_unused]] CSLConstList papszExtraOptions)
125	{
126	using namespace gdal;
127
128	viewshed::Options oOpts;	×
129	oOpts.outputFormat = pszDriverName;	×
130	oOpts.outputFilename = pszTargetRasterName;	×
131	oOpts.creationOpts = papszCreationOptions;	×
132	oOpts.observer.x = dfObserverX;	×
133	oOpts.observer.y = dfObserverY;	×
134	oOpts.observer.z = dfObserverHeight;	×
135	oOpts.targetHeight = dfTargetHeight;	×
136	oOpts.curveCoeff = dfCurvCoeff;	×
137	oOpts.maxDistance = dfMaxDistance;	×
138	oOpts.nodataVal = dfNoDataVal;	×
139
140	switch (eMode)	×
141	{
142	case GVM_Edge:	×
143	oOpts.cellMode = viewshed::CellMode::Edge;	×
144	break;	×
145	case GVM_Diagonal:	×
146	oOpts.cellMode = viewshed::CellMode::Diagonal;	×
147	break;	×
148	case GVM_Min:	×
149	oOpts.cellMode = viewshed::CellMode::Min;	×
150	break;	×
151	case GVM_Max:	×
152	oOpts.cellMode = viewshed::CellMode::Max;	×
153	break;	×
154	}
155
156	switch (heightMode)	×
157	{
158	case GVOT_MIN_TARGET_HEIGHT_FROM_DEM:	×
159	oOpts.outputMode = viewshed::OutputMode::DEM;	×
160	break;	×
161	case GVOT_MIN_TARGET_HEIGHT_FROM_GROUND:	×
162	oOpts.outputMode = viewshed::OutputMode::Ground;	×
163	break;	×
164	case GVOT_NORMAL:	×
165	oOpts.outputMode = viewshed::OutputMode::Normal;	×
166	break;	×
167	}
168
169	if (!GDALIsValueInRange<uint8_t>(dfVisibleVal))	×
170	{
171	CPLError(CE_Failure, CPLE_AppDefined,	×
172	"dfVisibleVal out of range. Must be [0, 255].");
173	return nullptr;	×
174	}
175	if (!GDALIsValueInRange<uint8_t>(dfInvisibleVal))	×
176	{
177	CPLError(CE_Failure, CPLE_AppDefined,	×
178	"dfInvisibleVal out of range. Must be [0, 255].");
179	return nullptr;	×
180	}
181	if (!GDALIsValueInRange<uint8_t>(dfOutOfRangeVal))	×
182	{
183	CPLError(CE_Failure, CPLE_AppDefined,	×
184	"dfOutOfRangeVal out of range. Must be [0, 255].");
185	return nullptr;	×
186	}
187	oOpts.visibleVal = dfVisibleVal;	×
188	oOpts.invisibleVal = dfInvisibleVal;	×
189	oOpts.outOfRangeVal = dfOutOfRangeVal;	×
190
191	gdal::viewshed::Viewshed v(oOpts);	×
192
193	if (!pfnProgress)	×
194	pfnProgress = GDALDummyProgress;	×
195	v.run(hBand, pfnProgress, pProgressArg);	×
196
197	return GDALDataset::FromHandle(v.output().release());	×
198	}
199
200	namespace gdal
201	{
202	namespace viewshed
203	{
204
205	namespace
206	{
207
208	bool getTransforms(GDALRasterBand &band, GDALGeoTransform &fwdTransform,	39✔
209	GDALGeoTransform &revTransform)
210	{
211	band.GetDataset()->GetGeoTransform(fwdTransform);	39✔
212	if (!fwdTransform.GetInverse(revTransform))	39✔
213	{
214	CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");	×
215	return false;	×
216	}
217	return true;	39✔
218	}
219
220	/// Shrink the extent of a window to just cover the slice defined by rays from
221	/// (nX, nY) and [startAngle, endAngle]
222	///
223	/// @param oOutExtent Window to modify
224	/// @param nX X coordinate of ray endpoint.
225	/// @param nY Y coordinate of ray endpoint.
226	/// @param startAngle Start angle of slice (standard mathmatics notion, in radians)
227	/// @param endAngle End angle of slice (standard mathmatics notion, in radians)
228	void shrinkWindowForAngles(Window &oOutExtent, int nX, int nY,	52✔
229	double startAngle, double endAngle)
230	{
231	/// NOTE: This probably doesn't work when the observer is outside the raster and
232	/// needs to be enhanced for that case.
233
234	if (startAngle == endAngle)	52✔
235	return;	37✔
236
237	Window win = oOutExtent;	15✔
238
239	// Set the X boundaries for the angles
240	int startAngleX = hIntersect(startAngle, nX, nY, win);	15✔
241	int stopAngleX = hIntersect(endAngle, nX, nY, win);	15✔
242
243	int xmax = nX;	15✔
244	if (!rayBetween(startAngle, endAngle, 0))	15✔
245	{
246	xmax = std::max(xmax, startAngleX);	7✔
247	xmax = std::max(xmax, stopAngleX);	7✔
248	// Add one to xmax since we want one past the stop. [start, stop)
249	oOutExtent.xStop = std::min(oOutExtent.xStop, xmax + 1);	7✔
250	}
251
252	int xmin = nX;	15✔
253	if (!rayBetween(startAngle, endAngle, M_PI))	15✔
254	{
255	xmin = std::min(xmin, startAngleX);	8✔
256	xmin = std::min(xmin, stopAngleX);	8✔
257	oOutExtent.xStart = std::max(oOutExtent.xStart, xmin);	8✔
258	}
259
260	// Set the Y boundaries for the angles
261	int startAngleY = vIntersect(startAngle, nX, nY, win);	15✔
262	int stopAngleY = vIntersect(endAngle, nX, nY, win);	15✔
263
264	int ymin = nY;	15✔
265	if (!rayBetween(startAngle, endAngle, M_PI / 2))	15✔
266	{
267	ymin = std::min(ymin, startAngleY);	7✔
268	ymin = std::min(ymin, stopAngleY);	7✔
269	oOutExtent.yStart = std::max(oOutExtent.yStart, ymin);	7✔
270	}
271	int ymax = nY;	15✔
272	if (!rayBetween(startAngle, endAngle, 3 * M_PI / 2))	15✔
273	{
274	ymax = std::max(ymax, startAngleY);	8✔
275	ymax = std::max(ymax, stopAngleY);	8✔
276	// Add one to ymax since we want one past the stop. [start, stop)
277	oOutExtent.yStop = std::min(oOutExtent.yStop, ymax + 1);	8✔
278	}
279	}
280
281	} // unnamed namespace
282
283	Viewshed::Viewshed(const Options &opts) : oOpts{opts}	39✔
284	{
285	}	39✔
286
287	Viewshed::~Viewshed() = default;
288
289	/// Calculate the extent of the output raster in terms of the input raster and
290	/// save the input raster extent.
291	///
292	/// @return false on error, true otherwise
293	bool Viewshed::calcExtents(int nX, int nY, const GDALGeoTransform &invGT)	39✔
294	{
295	// We start with the assumption that the output size matches the input.
296	oOutExtent.xStop = GDALGetRasterBandXSize(pSrcBand);	39✔
297	oOutExtent.yStop = GDALGetRasterBandYSize(pSrcBand);	39✔
298
299	if (!oOutExtent.contains(nX, nY))	39✔
300	{
301	if (oOpts.startAngle != oOpts.endAngle)	8✔
302	{
303	CPLError(CE_Failure, CPLE_AppDefined,	×
304	"Angle masking is not supported with an out-of-raster "
305	"observer.");
306	return false;	×
307	}
308	CPLError(CE_Warning, CPLE_AppDefined,	8✔
309	"NOTE: The observer location falls outside of the DEM area");
310	}
311
312	constexpr double EPSILON = 1e-8;	39✔
313	if (oOpts.maxDistance > 0)	39✔
314	{
315	//ABELL - This assumes that the transformation is only a scaling. Should be fixed.
316	// Find the distance in the direction of the transformed unit vector in the X and Y
317	// directions and use those factors to determine the limiting values in the raster space.
318	int nXStart = static_cast<int>(	3✔
319	std::floor(nX - invGT[1] * oOpts.maxDistance + EPSILON));	3✔
320	int nXStop = static_cast<int>(	3✔
321	std::ceil(nX + invGT[1] * oOpts.maxDistance - EPSILON) + 1);	3✔
322	//ABELL - These seem to be wrong. The transform of 1 is no transform, so not
323	// sure why we're adding one in the first case. Really, the transformed distance
324	// should add EPSILON. Not sure what the change should be for a negative transform,
325	// which is what I think is being handled with the 1/0 addition/subtraction.
326	int nYStart =
327	static_cast<int>(std::floor(	3✔
328	nY - std::fabs(invGT[5]) * oOpts.maxDistance + EPSILON)) -	3✔
329	(invGT[5] > 0 ? 1 : 0);	3✔
330	int nYStop = static_cast<int>(	3✔
331	std::ceil(nY + std::fabs(invGT[5]) * oOpts.maxDistance - EPSILON) +	6✔
332	(invGT[5] < 0 ? 1 : 0));	3✔
333
334	// If the limits are invalid, set the window size to zero to trigger the error below.
335	if (nXStart >= oOutExtent.xStop \|\| nXStop < 0 \|\|	3✔
336	nYStart >= oOutExtent.yStop \|\| nYStop < 0)	3✔
337	{
338	oOutExtent = Window();	×
339	}
340	else
341	{
342	oOutExtent.xStart = std::max(nXStart, 0);	3✔
343	oOutExtent.xStop = std::min(nXStop, oOutExtent.xStop);	3✔
344
345	oOutExtent.yStart = std::max(nYStart, 0);	3✔
346	oOutExtent.yStop = std::min(nYStop, oOutExtent.yStop);	3✔
347	}
348	}
349
350	if (oOutExtent.xSize() == 0 \|\| oOutExtent.ySize() == 0)	39✔
351	{
352	CPLError(CE_Failure, CPLE_AppDefined,	×
353	"Invalid target raster size due to transform "
354	"and/or distance limitation.");
355	return false;	×
356	}
357
358	shrinkWindowForAngles(oOutExtent, nX, nY, oOpts.startAngle, oOpts.endAngle);	39✔
359
360	// normalize horizontal index to [ 0, oOutExtent.xSize() )
361	oCurExtent = oOutExtent;	39✔
362	oCurExtent.shiftX(-oOutExtent.xStart);	39✔
363
364	return true;	39✔
365	}
366
367	/// Compute the viewshed of a raster band.
368	///
369	/// @param band Pointer to the raster band to be processed.
370	/// @param pfnProgress Pointer to the progress function. Can be null.
371	/// @param pProgressArg Argument passed to the progress function
372	/// @return True on success, false otherwise.
373	bool Viewshed::run(GDALRasterBandH band, GDALProgressFunc pfnProgress,	39✔
374	void *pProgressArg)
375	{
376	pSrcBand = static_cast<GDALRasterBand *>(band);	39✔
377
378	GDALGeoTransform fwdTransform, invTransform;	39✔
379	if (!getTransforms(*pSrcBand, fwdTransform, invTransform))	39✔
380	return false;	×
381
382	// calculate observer position
383	double dfX, dfY;
384	invTransform.Apply(oOpts.observer.x, oOpts.observer.y, &dfX, &dfY);	39✔
385	if (!GDALIsValueInRange<int>(dfX))	39✔
386	{
387	CPLError(CE_Failure, CPLE_AppDefined, "Observer X value out of range");	×
388	return false;	×
389	}
390	if (!GDALIsValueInRange<int>(dfY))	39✔
391	{
392	CPLError(CE_Failure, CPLE_AppDefined, "Observer Y value out of range");	×
393	return false;	×
394	}
395	int nX = static_cast<int>(dfX);	39✔
396	int nY = static_cast<int>(dfY);	39✔
397
398	if (oOpts.startAngle < 0 \|\| oOpts.startAngle >= 360)	39✔
399	{
400	CPLError(CE_Failure, CPLE_AppDefined,	×
401	"Start angle out of range. Must be [0, 360).");
402	return false;	×
403	}
404	if (oOpts.endAngle < 0 \|\| oOpts.endAngle >= 360)	39✔
405	{
406	CPLError(CE_Failure, CPLE_AppDefined,	×
407	"End angle out of range. Must be [0, 360).");
408	return false;	×
409	}
410	if (oOpts.highPitch > 90)	39✔
411	{
412	CPLError(CE_Failure, CPLE_AppDefined,	×
413	"Invalid highPitch. Cannot be greater than 90.");
414	return false;	×
415	}
416	if (oOpts.lowPitch < -90)	39✔
417	{
418	CPLError(CE_Failure, CPLE_AppDefined,	×
419	"Invalid lowPitch. Cannot be less than -90.");
420	return false;	×
421	}
422	if (oOpts.highPitch <= oOpts.lowPitch)	39✔
423	{
424	CPLError(CE_Failure, CPLE_AppDefined,	×
425	"Invalid pitch. highPitch must be > lowPitch");
426	return false;	×
427	}
428
429	// Normalize angle to radians and standard math arrangement.
430	oOpts.startAngle = normalizeAngle(oOpts.startAngle);	39✔
431	oOpts.endAngle = normalizeAngle(oOpts.endAngle);	39✔
432
433	// Must calculate extents in order to make the output dataset.
434	if (!calcExtents(nX, nY, invTransform))	39✔
435	return false;	×
436
437	poDstDS = createOutputDataset(*pSrcBand, oOpts, oOutExtent);	39✔
438	if (!poDstDS)	39✔
439	return false;	×
440
441	// Create the progress reporter.
442	Progress oProgress(pfnProgress, pProgressArg, oOutExtent.ySize());	78✔
443
444	// Execute the viewshed algorithm.
445	GDALRasterBand *pDstBand = poDstDS->GetRasterBand(1);	39✔
446	ViewshedExecutor executor(pSrcBand, pDstBand, nX, nY, oOutExtent,	39✔
447	oCurExtent, oOpts, oProgress,	39✔
448	/* emitWarningIfNoData = */ true);	78✔
449	executor.run();	39✔
450	oProgress.emit(1);	39✔
451	return static_cast<bool>(poDstDS);	39✔
452	}
453
454	// Adjust the coefficient of curvature for non-earth SRS.
455	/// \param curveCoeff Current curve coefficient
456	/// \param hSrcDS Source dataset
457	/// \return Adjusted curve coefficient.
458	double adjustCurveCoeff(double curveCoeff, GDALDatasetH hSrcDS)	14✔
459	{
460	const OGRSpatialReference *poSRS =
461	GDALDataset::FromHandle(hSrcDS)->GetSpatialRef();	14✔
462	if (poSRS)	14✔
463	{
464	OGRErr eSRSerr;
465	const double dfSemiMajor = poSRS->GetSemiMajor(&eSRSerr);	13✔
466	if (eSRSerr != OGRERR_FAILURE &&	13✔
467	fabs(dfSemiMajor - SRS_WGS84_SEMIMAJOR) >	13✔
468	0.05 * SRS_WGS84_SEMIMAJOR)
469	{
470	curveCoeff = 1.0;	×
471	CPLDebug("gdal_viewshed",	×
472	"Using -cc=1.0 as a non-Earth CRS has been detected");
473	}
474	}
475	return curveCoeff;	14✔
476	}
477
478	#if defined(__clang__) \|\| defined(__GNUC__)
479	#pragma GCC diagnostic ignored "-Wmissing-declarations"
480	#endif
481
482	void testShrinkWindowForAngles(Window &oOutExtent, int nX, int nY,	13✔
483	double startAngle, double endAngle)
484	{
485	shrinkWindowForAngles(oOutExtent, nX, nY, startAngle, endAngle);	13✔
486	}	13✔
487
488	} // namespace viewshed
489	} // namespace gdal

OSGeo / gdal / 15885686134

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous