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