Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Line of Sight
4 : * Purpose: Core algorithm implementation for line of sight algorithms.
5 : * Author: Ryan Friedman, ryanfriedman5410+gdal@gmail.com
6 : *
7 : ******************************************************************************
8 : *
9 : * SPDX-License-Identifier: MIT
10 : ****************************************************************************/
11 :
12 : #include <functional>
13 : #include <cmath>
14 :
15 : #include "cpl_port.h"
16 : #include "gdal_alg.h"
17 :
18 : // There's a plethora of bresenham implementations, all questionable production
19 : // quality. Bresenham optimizes for integer math, which makes sense for raster
20 : // datasets in 2D. For 3D, a 3D bresenham could be used if the altitude is also
21 : // integer resolution.
22 : // 2D:
23 : // https://codereview.stackexchange.com/questions/77460/bresenhams-line-algorithm-optimization
24 : // https://gist.github.com/ssavi-ict/092501c69e2ffec65e96a8865470ad2f
25 : // https://blog.demofox.org/2015/01/17/bresenhams-drawing-algorithms/
26 : // https://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm
27 : // https://www.cs.helsinki.fi/group/goa/mallinnus/lines/bresenh.html
28 : // https://stackoverflow.com/questions/10060046/drawing-lines-with-bresenhams-line-algorithm
29 : // http://www.edepot.com/linebresenham.html
30 : // 3D:
31 : // https://gist.github.com/yamamushi/5823518
32 :
33 : // Run bresenham terrain checking from (x1, y1) to (x2, y2).
34 : // The callback is run at every point along the line, which should return True
35 : // if the point is above terrain. Bresenham2D will return true if all points
36 : // have LOS between the start and end.
37 : static bool
38 13 : Bresenham2D(const int x1, const int y1, const int x2, const int y2,
39 : std::function<auto(const int, const int)->bool> OnBresenhamPoint)
40 : {
41 13 : bool isAboveTerrain = true;
42 : int dx, dy;
43 : int incx, incy;
44 :
45 13 : if (x2 >= x1)
46 : {
47 11 : dx = x2 - x1;
48 11 : incx = 1;
49 : }
50 : else
51 : {
52 2 : dx = x1 - x2;
53 2 : incx = -1;
54 : }
55 :
56 13 : if (y2 >= y1)
57 : {
58 8 : dy = y2 - y1;
59 8 : incy = 1;
60 : }
61 : else
62 : {
63 5 : dy = y1 - y2;
64 5 : incy = -1;
65 : }
66 :
67 13 : auto x = x1;
68 13 : auto y = y1;
69 : int balance;
70 :
71 13 : if (dx >= dy)
72 : {
73 10 : dy <<= 1;
74 10 : balance = dy - dx;
75 10 : dx *= 2;
76 :
77 290 : while (x != x2 && isAboveTerrain)
78 : {
79 280 : isAboveTerrain &= OnBresenhamPoint(x, y);
80 280 : if (balance >= 0)
81 : {
82 280 : y += incy;
83 280 : balance -= dx;
84 : }
85 280 : balance += dy;
86 280 : x += incx;
87 : }
88 10 : isAboveTerrain &= OnBresenhamPoint(x, y);
89 : }
90 : else
91 : {
92 3 : dx *= 2;
93 3 : balance = dx - dy;
94 3 : dy *= 2;
95 :
96 13 : while (y != y2 && isAboveTerrain)
97 : {
98 10 : isAboveTerrain &= OnBresenhamPoint(x, y);
99 10 : if (balance >= 0)
100 : {
101 6 : x += incx;
102 6 : balance -= dy;
103 : }
104 10 : balance += dx;
105 10 : y += incy;
106 : }
107 3 : isAboveTerrain &= OnBresenhamPoint(x, y);
108 : }
109 13 : return isAboveTerrain;
110 : }
111 :
112 : // Get the elevation of a single point.
113 421 : static bool GetElevation(const GDALRasterBandH hBand, const int x, const int y,
114 : double &val)
115 : {
116 : /// @todo GDALCachedPixelAccessor may give increased performance.
117 421 : return GDALRasterIO(hBand, GF_Read, x, y, 1, 1, &val, 1, 1, GDT_Float64, 0,
118 421 : 0) == CE_None;
119 : }
120 :
121 : // Check a single location is above (or equal) to terrain height.
122 421 : static bool IsAboveTerrain(const GDALRasterBandH hBand, const int x,
123 : const int y, const double z)
124 : {
125 : double terrainHeight;
126 421 : if (GetElevation(hBand, x, y, terrainHeight))
127 : {
128 420 : return z >= terrainHeight;
129 : }
130 : else
131 : {
132 1 : return false;
133 : }
134 : }
135 :
136 : /************************************************************************/
137 : /* GDALIsLineOfSightVisible() */
138 : /************************************************************************/
139 :
140 : /**
141 : * Check Line of Sight between two points.
142 : * Both input coordinates must be within the raster coordinate bounds.
143 : *
144 : * This algorithm will check line of sight using a Bresenham algorithm.
145 : * https://www.researchgate.net/publication/2411280_Efficient_Line-of-Sight_Algorithms_for_Real_Terrain_Data
146 : * Line of sight is computed in raster coordinate space, and thus may not be
147 : * appropriate. For example, datasets referenced against geographic coordinate
148 : * at high latitudes may have issues. A point exactly at the height of the DEM
149 : * is treated as visible from above.
150 : *
151 : * @param hBand The band to read the DEM data from. This must NOT be null.
152 : *
153 : * @param xA The X location (raster column) of the first point to check on the
154 : * raster.
155 : *
156 : * @param yA The Y location (raster row) of the first point to check on the
157 : * raster.
158 : *
159 : * @param zA The Z location (height) of the first point to check.
160 : *
161 : * @param xB The X location (raster column) of the second point to check on the
162 : * raster.
163 : *
164 : * @param yB The Y location (raster row) of the second point to check on the
165 : * raster.
166 : *
167 : * @param zB The Z location (height) of the second point to check.
168 : *
169 : * @param[out] pnxTerrainIntersection The X location where the LOS line
170 : * intersects with terrain, or nullptr if it does not intersect terrain.
171 : *
172 : * @param[out] pnyTerrainIntersection The Y location where the LOS line
173 : * intersects with terrain, or nullptr if it does not intersect terrain.
174 : *
175 : * @param papszOptions Options for the line of sight algorithm (currently
176 : * ignored).
177 : *
178 : * @return True if the two points are within Line of Sight.
179 : *
180 : * @since GDAL 3.9
181 : */
182 :
183 37 : bool GDALIsLineOfSightVisible(const GDALRasterBandH hBand, const int xA,
184 : const int yA, const double zA, const int xB,
185 : const int yB, const double zB,
186 : int *pnxTerrainIntersection,
187 : int *pnyTerrainIntersection,
188 : CPL_UNUSED CSLConstList papszOptions)
189 : {
190 37 : VALIDATE_POINTER1(hBand, "GDALIsLineOfSightVisible", false);
191 :
192 : // A lambda to set the X-Y intersection if it's not null
193 28 : auto SetXYIntersection = [&](const int x, const int y)
194 : {
195 28 : if (pnxTerrainIntersection != nullptr)
196 : {
197 10 : *pnxTerrainIntersection = x;
198 : }
199 28 : if (pnyTerrainIntersection != nullptr)
200 : {
201 10 : *pnyTerrainIntersection = y;
202 : }
203 65 : };
204 :
205 37 : if (pnxTerrainIntersection)
206 12 : *pnxTerrainIntersection = -1;
207 :
208 37 : if (pnyTerrainIntersection)
209 12 : *pnyTerrainIntersection = -1;
210 :
211 : // Perform a preliminary check of the start and end points.
212 37 : if (!IsAboveTerrain(hBand, xA, yA, zA))
213 : {
214 10 : SetXYIntersection(xA, yA);
215 10 : return false;
216 : }
217 27 : if (!IsAboveTerrain(hBand, xB, yB, zB))
218 : {
219 2 : SetXYIntersection(xB, yB);
220 2 : return false;
221 : }
222 :
223 : // If both X and Y are the same, no further checks are needed.
224 25 : if (xA == xB && yA == yB)
225 : {
226 3 : return true;
227 : }
228 :
229 : // Lambda for Linear interpolate like C++20 std::lerp.
230 357 : auto lerp = [](const double a, const double b, const double t)
231 357 : { return a + t * (b - a); };
232 :
233 : // Lambda for getting Z test height given y input along the LOS line.
234 : // Only to be used for vertical line checks.
235 18 : auto GetZValueFromY = [&](const int y) -> double
236 : {
237 : // A ratio of 0.0 corresponds to being at yA.
238 18 : const auto ratio =
239 18 : static_cast<double>(y - yA) / static_cast<double>(yB - yA);
240 18 : return lerp(zA, zB, ratio);
241 22 : };
242 :
243 : // Lambda for getting Z test height given x input along the LOS line.
244 : // Only to be used for horizontal line checks.
245 36 : auto GetZValueFromX = [&](const int x) -> double
246 : {
247 : // A ratio of 0.0 corresponds to being at xA.
248 36 : const auto ratio =
249 36 : static_cast<double>(x - xA) / static_cast<double>(xB - xA);
250 36 : return lerp(zA, zB, ratio);
251 22 : };
252 :
253 : // Lambda for checking path safety of a vertical line.
254 : // Returns true if the path has clear LOS.
255 4 : auto CheckVerticalLine = [&]() -> bool
256 : {
257 4 : CPLAssert(xA == xB);
258 4 : CPLAssert(yA != yB);
259 :
260 4 : if (yA < yB)
261 : {
262 10 : for (int y = yA; y <= yB; ++y)
263 : {
264 9 : const auto zTest = GetZValueFromY(y);
265 9 : if (!IsAboveTerrain(hBand, xA, y, zTest))
266 : {
267 1 : SetXYIntersection(xA, y);
268 1 : return false;
269 : }
270 : }
271 1 : return true;
272 : }
273 : else
274 : {
275 10 : for (int y = yA; y >= yB; --y)
276 : {
277 9 : const auto zTest = GetZValueFromY(y);
278 9 : if (!IsAboveTerrain(hBand, xA, y, zTest))
279 : {
280 1 : SetXYIntersection(xA, y);
281 1 : return false;
282 : }
283 : }
284 1 : return true;
285 : }
286 22 : };
287 :
288 : // Lambda for checking path safety of a horizontal line.
289 : // Returns true if the path has clear LOS.
290 5 : auto CheckHorizontalLine = [&]() -> bool
291 : {
292 5 : CPLAssert(yA == yB);
293 5 : CPLAssert(xA != xB);
294 :
295 5 : if (xA < xB)
296 : {
297 21 : for (int x = xA; x <= xB; ++x)
298 : {
299 19 : const auto zTest = GetZValueFromX(x);
300 19 : if (!IsAboveTerrain(hBand, x, yA, zTest))
301 : {
302 1 : SetXYIntersection(x, yA);
303 1 : return false;
304 : }
305 : }
306 2 : return true;
307 : }
308 : else
309 : {
310 18 : for (int x = xA; x >= xB; --x)
311 : {
312 17 : const auto zTest = GetZValueFromX(x);
313 17 : if (!IsAboveTerrain(hBand, x, yA, zTest))
314 : {
315 1 : SetXYIntersection(x, yA);
316 1 : return false;
317 : }
318 : }
319 1 : return true;
320 : }
321 22 : };
322 :
323 : // Handle special cases if it's a vertical or horizontal line
324 : // (don't use bresenham).
325 22 : if (xA == xB)
326 : {
327 4 : return CheckVerticalLine();
328 : }
329 18 : if (yA == yB)
330 : {
331 5 : return CheckHorizontalLine();
332 : }
333 :
334 : // Use an interpolated Z height with 2D bresenham for the remaining cases.
335 :
336 : // Lambda for computing the square of a number
337 1212 : auto SQUARE = [](const double d) -> double { return d * d; };
338 :
339 : // Lambda for getting Z test height given x-y input along bresenham line.
340 303 : auto GetZValueFromXY = [&](const int x, const int y) -> double
341 : {
342 303 : const auto rNum = SQUARE(static_cast<double>(x - xA)) +
343 303 : SQUARE(static_cast<double>(y - yA));
344 303 : const auto rDenom = SQUARE(static_cast<double>(xB - xA)) +
345 303 : SQUARE(static_cast<double>(yB - yA));
346 : /// @todo In order to reduce CPU cost and avoid a sqrt operation,
347 : /// consider the approach to just the ratio along x or y depending on
348 : /// whether the line is steep or shallow.
349 : /// See https://github.com/OSGeo/gdal/pull/9506#discussion_r1532459689.
350 : const double ratio =
351 303 : sqrt(static_cast<double>(rNum) / static_cast<double>(rDenom));
352 303 : return lerp(zA, zB, ratio);
353 13 : };
354 :
355 : // Lambda to get if above terrain at a bresenham-computed location.
356 303 : auto OnBresenhamPoint = [&](const int x, const int y) -> bool
357 : {
358 303 : const auto z = GetZValueFromXY(x, y);
359 303 : const auto isAbove = IsAboveTerrain(hBand, x, y, z);
360 303 : if (!isAbove)
361 : {
362 12 : SetXYIntersection(x, y);
363 : }
364 303 : return isAbove;
365 13 : };
366 :
367 13 : return Bresenham2D(xA, yA, xB, yB, OnBresenhamPoint);
368 : }
|
