/src/gdal/alg/viewshed/util.cpp

Source
/******************************************************************************
 * (c) 2024 info@hobu.co
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include <array>
#include <cmath>
#include <iostream>
#include <limits>

#include "gdal_priv.h"
#include "util.h"
#include "viewshed_types.h"

namespace gdal
{
namespace viewshed
{

/// Normalize a masking angle. Change from clockwise with 0 north (up) to counterclockwise
/// with 0 to the east (right) and change to radians.
///
// @param maskAngle  Masking angle in degrees.
double normalizeAngle(double maskAngle)
{
    maskAngle = 90 - maskAngle;
    if (maskAngle < 0)
        maskAngle += 360;
    return maskAngle * (M_PI / 180);
}

/// Compute the X intersect position on the line Y = y with a ray extending
/// from (nX, nY) along `angle`.
///
/// @param angle  Angle in radians, standard arrangement.
/// @param nX  X coordinate of ray endpoint.
/// @param nY  Y coordinate of ray endpoint.
/// @param y  Horizontal line where Y = y.
/// @return  X intersect or NaN
double horizontalIntersect(double angle, int nX, int nY, int y)
{
    double x = std::numeric_limits<double>::quiet_NaN();

    if (nY == y)
        x = nX;
    else if (nY > y)
    {
        if (ARE_REAL_EQUAL(angle, M_PI / 2))
            x = nX;
        else if (angle > 0 && angle < M_PI)
            x = nX + (nY - y) / std::tan(angle);
    }
    else  // nY < y
    {
        if (ARE_REAL_EQUAL(angle, 3 * M_PI / 2))
            x = nX;
        else if (angle > M_PI)
            x = nX - (y - nY) / std::tan(angle);
    }
    return x;
}

/// Compute the X intersect position on the line Y = y with a ray extending
/// from (nX, nY) along `angle`.
///
/// @param angle  Angle in radians, standard arrangement.
/// @param nX  X coordinate of ray endpoint.
/// @param nY  Y coordinate of ray endpoint.
/// @param y  Horizontal line where Y = y.
/// @return  Rounded X intersection of the sentinel INVALID_ISECT
int hIntersect(double angle, int nX, int nY, int y)
{
    double x = horizontalIntersect(angle, nX, nY, y);
    if (std::isnan(x))
        return INVALID_ISECT;
    return static_cast<int>(std::round(x));
}

/// Compute the X intersect on one of the horizontal edges of a window
/// with a ray extending from (nX, nY) along `angle`, clamped the the extent of a window.
///
/// @param angle  Angle in radians, standard arrangement.
/// @param nX  X coordinate of ray endpoint.
/// @param nY  Y coordinate of ray endpoint.
/// @param win  Window to intersect.
/// @return  X intersect, clamped to the window extent.
int hIntersect(double angle, int nX, int nY, const Window &win)
{
    if (ARE_REAL_EQUAL(angle, M_PI))
        return win.xStart;
    if (ARE_REAL_EQUAL(angle, 0))
        return win.xStop;
    double x = horizontalIntersect(angle, nX, nY, win.yStart);
    if (std::isnan(x))
        x = horizontalIntersect(angle, nX, nY, win.yStop);
    return std::clamp(static_cast<int>(std::round(x)), win.xStart, win.xStop);
}

/// Compute the X intersect position on the line Y = y with a ray extending
/// from (nX, nY) along `angle`.
///
/// @param angle  Angle in radians, standard arrangement.
/// @param nX  X coordinate of ray endpoint.
/// @param nY  Y coordinate of ray endpoint.
/// @param x  Vertical  line where X = x.
/// @return  Y intersect or NaN
double verticalIntersect(double angle, int nX, int nY, int x)
{
    double y = std::numeric_limits<double>::quiet_NaN();

    if (nX == x)
        y = nY;
    else if (nX < x)
    {
        if (ARE_REAL_EQUAL(angle, 0))
            y = nY;
        else if (angle < M_PI / 2 || angle > M_PI * 3 / 2)
            y = nY + (nX - x) * std::tan(angle);
    }
    else  // nX > x
    {
        if (ARE_REAL_EQUAL(angle, M_PI))
            y = nY;
        else if (angle > M_PI / 2 && angle < M_PI * 3 / 2)
            y = nY - (x - nX) * std::tan(angle);
    }
    return y;
}

/// Compute the X intersect position on the line Y = y with a ray extending
/// from (nX, nY) along `angle`.
///
/// @param angle  Angle in radians, standard arrangement.
/// @param nX  X coordinate of ray endpoint.
/// @param nY  Y coordinate of ray endpoint.
/// @param x  Horizontal line where X = x.
/// @return  Rounded Y intersection of the sentinel INVALID_ISECT
int vIntersect(double angle, int nX, int nY, int x)
{
    double y = verticalIntersect(angle, nX, nY, x);
    if (std::isnan(y))
        return INVALID_ISECT;
    return static_cast<int>(std::round(y));
}

/// Compute the Y intersect on one of the vertical edges of a window
/// with a ray extending from (nX, nY) along `angle`, clamped the the extent
/// of the window.
///
/// @param angle  Angle in radians, standard arrangement.
/// @param nX  X coordinate of ray endpoint.
/// @param nY  Y coordinate of ray endpoint.
/// @param win  Window to intersect.
/// @return  y intersect, clamped to the window extent.
int vIntersect(double angle, int nX, int nY, const Window &win)
{
    if (ARE_REAL_EQUAL(angle, M_PI / 2))
        return win.yStart;
    if (ARE_REAL_EQUAL(angle, 3 * M_PI / 2))
        return win.yStop;
    double y = verticalIntersect(angle, nX, nY, win.xStart);
    if (std::isnan(y))
        y = verticalIntersect(angle, nX, nY, win.xStop);
    return std::clamp(static_cast<int>(std::round(y)), win.yStart, win.yStop);
}

/// Determine if ray is in the slice between two rays starting at `start` and
/// going clockwise to `end` (inclusive). [start, end]
/// @param  start  Start angle.
/// @param  end  End angle.
/// @param  test  Test angle.
/// @return  Whether `test` lies in the slice [start, end]
bool rayBetween(double start, double end, double test)
{
    // Our angles go counterclockwise, so swap start and end
    std::swap(start, end);
    if (start < end)
        return (test >= start && test <= end);
    else if (start > end)
        return (test >= start || test <= end);
    return false;
}

/// Get the band size
///
/// @param  band Raster band
/// @return  The raster band size.
size_t bandSize(GDALRasterBand &band)
{
    return static_cast<size_t>(band.GetXSize()) * band.GetYSize();
}

/// Create the output dataset.
///
/// @param  srcBand  Source raster band.
/// @param  opts  Options.
/// @param  extent  Output dataset extent.
/// @return  The output dataset to be filled with data.
DatasetPtr createOutputDataset(GDALRasterBand &srcBand, const Options &opts,
                               const Window &extent)
{
    GDALDriverManager *hMgr = GetGDALDriverManager();
    GDALDriver *hDriver = hMgr->GetDriverByName(opts.outputFormat.c_str());
    if (!hDriver)
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Cannot get driver");
        return nullptr;
    }

    /* create output raster */
    DatasetPtr dataset(hDriver->Create(
        opts.outputFilename.c_str(), extent.xSize(), extent.ySize(), 1,
        opts.outputMode == OutputMode::Normal ? GDT_UInt8 : GDT_Float64,
        const_cast<char **>(opts.creationOpts.List())));
    if (!dataset)
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Cannot create dataset for %s",
                 opts.outputFilename.c_str());
        return nullptr;
    }

    /* copy srs */
    dataset->SetSpatialRef(srcBand.GetDataset()->GetSpatialRef());

    GDALGeoTransform srcGT, dstGT;
    srcBand.GetDataset()->GetGeoTransform(srcGT);
    dstGT[0] = srcGT[0] + srcGT[1] * extent.xStart + srcGT[2] * extent.yStart;
    dstGT[1] = srcGT[1];
    dstGT[2] = srcGT[2];
    dstGT[3] = srcGT[3] + srcGT[4] * extent.xStart + srcGT[5] * extent.yStart;
    dstGT[4] = srcGT[4];
    dstGT[5] = srcGT[5];
    dataset->SetGeoTransform(dstGT);

    GDALRasterBand *pBand = dataset->GetRasterBand(1);
    if (!pBand)
    {
        CPLError(CE_Failure, CPLE_AppDefined, "Cannot get band for %s",
                 opts.outputFilename.c_str());
        return nullptr;
    }

    if (opts.nodataVal >= 0)
        GDALSetRasterNoDataValue(pBand, opts.nodataVal);
    return dataset;
}

}  // namespace viewshed
}  // namespace gdal

Coverage Report

Created: 2026-02-14 09:00

Line	Count	Source
1		/******************************************************************************
2		* (c) 2024 info@hobu.co
3		*
4		* SPDX-License-Identifier: MIT
5		****************************************************************************/
6
7		#include <array>
8		#include <cmath>
9		#include <iostream>
10		#include <limits>
11
12		#include "gdal_priv.h"
13		#include "util.h"
14		#include "viewshed_types.h"
15
16		namespace gdal
17		{
18		namespace viewshed
19		{
20
21		/// Normalize a masking angle. Change from clockwise with 0 north (up) to counterclockwise
22		/// with 0 to the east (right) and change to radians.
23		///
24		// @param maskAngle Masking angle in degrees.
25		double normalizeAngle(double maskAngle)
26	0	{
27	0	maskAngle = 90 - maskAngle;
28	0	if (maskAngle < 0)
29	0	maskAngle += 360;
30	0	return maskAngle * (M_PI / 180);
31	0	}
32
33		/// Compute the X intersect position on the line Y = y with a ray extending
34		/// from (nX, nY) along `angle`.
35		///
36		/// @param angle Angle in radians, standard arrangement.
37		/// @param nX X coordinate of ray endpoint.
38		/// @param nY Y coordinate of ray endpoint.
39		/// @param y Horizontal line where Y = y.
40		/// @return X intersect or NaN
41		double horizontalIntersect(double angle, int nX, int nY, int y)
42	0	{
43	0	double x = std::numeric_limits<double>::quiet_NaN();
44
45	0	if (nY == y)
46	0	x = nX;
47	0	else if (nY > y)
48	0	{
49	0	if (ARE_REAL_EQUAL(angle, M_PI / 2))
50	0	x = nX;
51	0	else if (angle > 0 && angle < M_PI)
52	0	x = nX + (nY - y) / std::tan(angle);
53	0	}
54	0	else // nY < y
55	0	{
56	0	if (ARE_REAL_EQUAL(angle, 3 * M_PI / 2))
57	0	x = nX;
58	0	else if (angle > M_PI)
59	0	x = nX - (y - nY) / std::tan(angle);
60	0	}
61	0	return x;
62	0	}
63
64		/// Compute the X intersect position on the line Y = y with a ray extending
65		/// from (nX, nY) along `angle`.
66		///
67		/// @param angle Angle in radians, standard arrangement.
68		/// @param nX X coordinate of ray endpoint.
69		/// @param nY Y coordinate of ray endpoint.
70		/// @param y Horizontal line where Y = y.
71		/// @return Rounded X intersection of the sentinel INVALID_ISECT
72		int hIntersect(double angle, int nX, int nY, int y)
73	0	{
74	0	double x = horizontalIntersect(angle, nX, nY, y);
75	0	if (std::isnan(x))
76	0	return INVALID_ISECT;
77	0	return static_cast<int>(std::round(x));
78	0	}
79
80		/// Compute the X intersect on one of the horizontal edges of a window
81		/// with a ray extending from (nX, nY) along `angle`, clamped the the extent of a window.
82		///
83		/// @param angle Angle in radians, standard arrangement.
84		/// @param nX X coordinate of ray endpoint.
85		/// @param nY Y coordinate of ray endpoint.
86		/// @param win Window to intersect.
87		/// @return X intersect, clamped to the window extent.
88		int hIntersect(double angle, int nX, int nY, const Window &win)
89	0	{
90	0	if (ARE_REAL_EQUAL(angle, M_PI))
91	0	return win.xStart;
92	0	if (ARE_REAL_EQUAL(angle, 0))
93	0	return win.xStop;
94	0	double x = horizontalIntersect(angle, nX, nY, win.yStart);
95	0	if (std::isnan(x))
96	0	x = horizontalIntersect(angle, nX, nY, win.yStop);
97	0	return std::clamp(static_cast<int>(std::round(x)), win.xStart, win.xStop);
98	0	}
99
100		/// Compute the X intersect position on the line Y = y with a ray extending
101		/// from (nX, nY) along `angle`.
102		///
103		/// @param angle Angle in radians, standard arrangement.
104		/// @param nX X coordinate of ray endpoint.
105		/// @param nY Y coordinate of ray endpoint.
106		/// @param x Vertical line where X = x.
107		/// @return Y intersect or NaN
108		double verticalIntersect(double angle, int nX, int nY, int x)
109	0	{
110	0	double y = std::numeric_limits<double>::quiet_NaN();
111
112	0	if (nX == x)
113	0	y = nY;
114	0	else if (nX < x)
115	0	{
116	0	if (ARE_REAL_EQUAL(angle, 0))
117	0	y = nY;
118	0	else if (angle < M_PI / 2 \|\| angle > M_PI * 3 / 2)
119	0	y = nY + (nX - x) * std::tan(angle);
120	0	}
121	0	else // nX > x
122	0	{
123	0	if (ARE_REAL_EQUAL(angle, M_PI))
124	0	y = nY;
125	0	else if (angle > M_PI / 2 && angle < M_PI * 3 / 2)
126	0	y = nY - (x - nX) * std::tan(angle);
127	0	}
128	0	return y;
129	0	}
130
131		/// Compute the X intersect position on the line Y = y with a ray extending
132		/// from (nX, nY) along `angle`.
133		///
134		/// @param angle Angle in radians, standard arrangement.
135		/// @param nX X coordinate of ray endpoint.
136		/// @param nY Y coordinate of ray endpoint.
137		/// @param x Horizontal line where X = x.
138		/// @return Rounded Y intersection of the sentinel INVALID_ISECT
139		int vIntersect(double angle, int nX, int nY, int x)
140	0	{
141	0	double y = verticalIntersect(angle, nX, nY, x);
142	0	if (std::isnan(y))
143	0	return INVALID_ISECT;
144	0	return static_cast<int>(std::round(y));
145	0	}
146
147		/// Compute the Y intersect on one of the vertical edges of a window
148		/// with a ray extending from (nX, nY) along `angle`, clamped the the extent
149		/// of the window.
150		///
151		/// @param angle Angle in radians, standard arrangement.
152		/// @param nX X coordinate of ray endpoint.
153		/// @param nY Y coordinate of ray endpoint.
154		/// @param win Window to intersect.
155		/// @return y intersect, clamped to the window extent.
156		int vIntersect(double angle, int nX, int nY, const Window &win)
157	0	{
158	0	if (ARE_REAL_EQUAL(angle, M_PI / 2))
159	0	return win.yStart;
160	0	if (ARE_REAL_EQUAL(angle, 3 * M_PI / 2))
161	0	return win.yStop;
162	0	double y = verticalIntersect(angle, nX, nY, win.xStart);
163	0	if (std::isnan(y))
164	0	y = verticalIntersect(angle, nX, nY, win.xStop);
165	0	return std::clamp(static_cast<int>(std::round(y)), win.yStart, win.yStop);
166	0	}
167
168		/// Determine if ray is in the slice between two rays starting at `start` and
169		/// going clockwise to `end` (inclusive). [start, end]
170		/// @param start Start angle.
171		/// @param end End angle.
172		/// @param test Test angle.
173		/// @return Whether `test` lies in the slice [start, end]
174		bool rayBetween(double start, double end, double test)
175	0	{
176		// Our angles go counterclockwise, so swap start and end
177	0	std::swap(start, end);
178	0	if (start < end)
179	0	return (test >= start && test <= end);
180	0	else if (start > end)
181	0	return (test >= start \|\| test <= end);
182	0	return false;
183	0	}
184
185		/// Get the band size
186		///
187		/// @param band Raster band
188		/// @return The raster band size.
189		size_t bandSize(GDALRasterBand &band)
190	0	{
191	0	return static_cast<size_t>(band.GetXSize()) * band.GetYSize();
192	0	}
193
194		/// Create the output dataset.
195		///
196		/// @param srcBand Source raster band.
197		/// @param opts Options.
198		/// @param extent Output dataset extent.
199		/// @return The output dataset to be filled with data.
200		DatasetPtr createOutputDataset(GDALRasterBand &srcBand, const Options &opts,
201		const Window &extent)
202	0	{
203	0	GDALDriverManager *hMgr = GetGDALDriverManager();
204	0	GDALDriver *hDriver = hMgr->GetDriverByName(opts.outputFormat.c_str());
205	0	if (!hDriver)
206	0	{
207	0	CPLError(CE_Failure, CPLE_AppDefined, "Cannot get driver");
208	0	return nullptr;
209	0	}
210
211		/* create output raster */
212	0	DatasetPtr dataset(hDriver->Create(
213	0	opts.outputFilename.c_str(), extent.xSize(), extent.ySize(), 1,
214	0	opts.outputMode == OutputMode::Normal ? GDT_UInt8 : GDT_Float64,
215	0	const_cast<char **>(opts.creationOpts.List())));
216	0	if (!dataset)
217	0	{
218	0	CPLError(CE_Failure, CPLE_AppDefined, "Cannot create dataset for %s",
219	0	opts.outputFilename.c_str());
220	0	return nullptr;
221	0	}
222
223		/* copy srs */
224	0	dataset->SetSpatialRef(srcBand.GetDataset()->GetSpatialRef());
225
226	0	GDALGeoTransform srcGT, dstGT;
227	0	srcBand.GetDataset()->GetGeoTransform(srcGT);
228	0	dstGT[0] = srcGT[0] + srcGT[1] * extent.xStart + srcGT[2] * extent.yStart;
229	0	dstGT[1] = srcGT[1];
230	0	dstGT[2] = srcGT[2];
231	0	dstGT[3] = srcGT[3] + srcGT[4] * extent.xStart + srcGT[5] * extent.yStart;
232	0	dstGT[4] = srcGT[4];
233	0	dstGT[5] = srcGT[5];
234	0	dataset->SetGeoTransform(dstGT);
235
236	0	GDALRasterBand *pBand = dataset->GetRasterBand(1);
237	0	if (!pBand)
238	0	{
239	0	CPLError(CE_Failure, CPLE_AppDefined, "Cannot get band for %s",
240	0	opts.outputFilename.c_str());
241	0	return nullptr;
242	0	}
243
244	0	if (opts.nodataVal >= 0)
245	0	GDALSetRasterNoDataValue(pBand, opts.nodataVal);
246	0	return dataset;
247	0	}
248
249		} // namespace viewshed
250		} // namespace gdal