Line data Source code
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 84 : double normalizeAngle(double maskAngle)
26 : {
27 84 : maskAngle = 90 - maskAngle;
28 84 : if (maskAngle < 0)
29 4 : maskAngle += 360;
30 84 : return maskAngle * (M_PI / 180);
31 : }
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 113 : double horizontalIntersect(double angle, int nX, int nY, int y)
42 : {
43 113 : double x = std::numeric_limits<double>::quiet_NaN();
44 :
45 113 : if (nY == y)
46 0 : x = nX;
47 113 : else if (nY > y)
48 : {
49 74 : if (ARE_REAL_EQUAL(angle, M_PI / 2))
50 21 : x = nX;
51 53 : else if (angle > 0 && angle < M_PI)
52 38 : x = nX + (nY - y) / std::tan(angle);
53 : }
54 : else // nY < y
55 : {
56 39 : if (ARE_REAL_EQUAL(angle, 3 * M_PI / 2))
57 1 : x = nX;
58 38 : else if (angle > M_PI)
59 14 : x = nX - (y - nY) / std::tan(angle);
60 : }
61 113 : return x;
62 : }
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 56 : int hIntersect(double angle, int nX, int nY, int y)
73 : {
74 56 : double x = horizontalIntersect(angle, nX, nY, y);
75 56 : if (std::isnan(x))
76 20 : return INVALID_ISECT;
77 36 : return static_cast<int>(std::round(x));
78 : }
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 30 : int hIntersect(double angle, int nX, int nY, const Window &win)
89 : {
90 30 : if (ARE_REAL_EQUAL(angle, M_PI))
91 0 : return win.xStart;
92 30 : if (ARE_REAL_EQUAL(angle, 0))
93 0 : return win.xStop;
94 30 : double x = horizontalIntersect(angle, nX, nY, win.yStart);
95 30 : if (std::isnan(x))
96 11 : x = horizontalIntersect(angle, nX, nY, win.yStop);
97 30 : return std::clamp(static_cast<int>(std::round(x)), win.xStart, win.xStop);
98 : }
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 61 : double verticalIntersect(double angle, int nX, int nY, int x)
109 : {
110 61 : double y = std::numeric_limits<double>::quiet_NaN();
111 :
112 61 : if (nX == x)
113 0 : y = nY;
114 61 : else if (nX < x)
115 : {
116 26 : if (ARE_REAL_EQUAL(angle, 0))
117 1 : y = nY;
118 25 : else if (angle < M_PI / 2 || angle > M_PI * 3 / 2)
119 20 : y = nY + (nX - x) * std::tan(angle);
120 : }
121 : else // nX > x
122 : {
123 35 : if (ARE_REAL_EQUAL(angle, M_PI))
124 1 : y = nY;
125 34 : else if (angle > M_PI / 2 && angle < M_PI * 3 / 2)
126 14 : y = nY - (x - nX) * std::tan(angle);
127 : }
128 61 : return y;
129 : }
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 0 : int vIntersect(double angle, int nX, int nY, int x)
140 : {
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 : }
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 30 : int vIntersect(double angle, int nX, int nY, const Window &win)
157 : {
158 30 : if (ARE_REAL_EQUAL(angle, M_PI / 2))
159 2 : return win.yStart;
160 28 : if (ARE_REAL_EQUAL(angle, 3 * M_PI / 2))
161 0 : return win.yStop;
162 28 : double y = verticalIntersect(angle, nX, nY, win.xStart);
163 28 : if (std::isnan(y))
164 17 : y = verticalIntersect(angle, nX, nY, win.xStop);
165 28 : return std::clamp(static_cast<int>(std::round(y)), win.yStart, win.yStop);
166 : }
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 103 : bool rayBetween(double start, double end, double test)
175 : {
176 : // Our angles go counterclockwise, so swap start and end
177 103 : std::swap(start, end);
178 103 : if (start < end)
179 39 : return (test >= start && test <= end);
180 64 : else if (start > end)
181 64 : return (test >= start || test <= end);
182 0 : return false;
183 : }
184 :
185 : /// Get the band size
186 : ///
187 : /// @param band Raster band
188 : /// @return The raster band size.
189 190 : size_t bandSize(GDALRasterBand &band)
190 : {
191 190 : return static_cast<size_t>(band.GetXSize()) * band.GetYSize();
192 : }
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 41 : DatasetPtr createOutputDataset(GDALRasterBand &srcBand, const Options &opts,
201 : const Window &extent)
202 : {
203 41 : GDALDriverManager *hMgr = GetGDALDriverManager();
204 41 : GDALDriver *hDriver = hMgr->GetDriverByName(opts.outputFormat.c_str());
205 41 : if (!hDriver)
206 : {
207 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot get driver");
208 0 : return nullptr;
209 : }
210 :
211 : /* create output raster */
212 : DatasetPtr dataset(hDriver->Create(
213 : opts.outputFilename.c_str(), extent.xSize(), extent.ySize(), 1,
214 41 : opts.outputMode == OutputMode::Normal ? GDT_Byte : GDT_Float64,
215 82 : const_cast<char **>(opts.creationOpts.List())));
216 41 : if (!dataset)
217 : {
218 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot create dataset for %s",
219 : opts.outputFilename.c_str());
220 0 : return nullptr;
221 : }
222 :
223 : /* copy srs */
224 41 : dataset->SetSpatialRef(srcBand.GetDataset()->GetSpatialRef());
225 :
226 : std::array<double, 6> adfSrcTransform;
227 : std::array<double, 6> adfDstTransform;
228 41 : srcBand.GetDataset()->GetGeoTransform(adfSrcTransform.data());
229 41 : adfDstTransform[0] = adfSrcTransform[0] +
230 41 : adfSrcTransform[1] * extent.xStart +
231 41 : adfSrcTransform[2] * extent.yStart;
232 41 : adfDstTransform[1] = adfSrcTransform[1];
233 41 : adfDstTransform[2] = adfSrcTransform[2];
234 41 : adfDstTransform[3] = adfSrcTransform[3] +
235 41 : adfSrcTransform[4] * extent.xStart +
236 41 : adfSrcTransform[5] * extent.yStart;
237 41 : adfDstTransform[4] = adfSrcTransform[4];
238 41 : adfDstTransform[5] = adfSrcTransform[5];
239 41 : dataset->SetGeoTransform(adfDstTransform.data());
240 :
241 41 : GDALRasterBand *pBand = dataset->GetRasterBand(1);
242 41 : if (!pBand)
243 : {
244 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot get band for %s",
245 : opts.outputFilename.c_str());
246 0 : return nullptr;
247 : }
248 :
249 41 : if (opts.nodataVal >= 0)
250 1 : GDALSetRasterNoDataValue(pBand, opts.nodataVal);
251 41 : return dataset;
252 : }
253 :
254 : } // namespace viewshed
255 : } // namespace gdal
|
