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