speeding up r.horizon

[girishnand]

@petrasovaa
(Please tag anyone else who may find this of interest.)

r.horizon, in its current incarnation uses an algorithm which has O(n^3) time complexity in the worst case, where n is the number of points along the side of a square raster. It does somewhat better than O(n^3) in typical use cases by using the following two tricks -

1. Maximum height is computed for each 100x100 cell block and the entire block is passed over if this height is less than the minimal height, over that block, of a ray from a point to its current estimated horizon.
2. If the height of a ray from a line towards its current estimated horizon exceeds the maximum pixel height in the raster, the remaining pixels in that direction are not checked.

But these tricks won't help much if you are computing the horizon towards a gradually rising slope, for example.

The following purportedly more efficient algorithms are available in the literature -

1. The Dozier and Bruno algorithm of 1981 has O(n^2) time complexity for each horizon direction and appears to be supremely simpler to implement as per code snippets included in the following publications -

• A Faster Solution to the Horizon Problem, J Dozier & J Bruno (1981), Comput & Geosci

• Rapid Calculation of Terrain Parameters for Radiation Modeling from Digital Elevation Data, Jeff Dozier & James Frew (1990), IEEE Trans on Geosci & Remote Sens

• Revisiting Topographic Horizons in the Era of Big Data and Parallel Computing, Jeff Dozier (2022), IEEE Geosci & Remote Sensing Lettrs

  It also appears to be widely used and opensource implementations appear to be available at the following locations -
  https://in.mathworks.com/matlabcentral/fileexchange/94800-topographic-horizons
  https://github.com/USDA-ARS-NWRC/topocalc/tree/main
  https://github.com/maximlamare/REDRESS/blob/main/redress/topography/horizon.py

2. The 1997 algorithm of Stewart purports to have O(n^2 (log(n))^2 ) time complexity for each horizon direction and is described in the following publication -

• Fast Horizon Computation at All Points of a Terrain with Visibility and Shading Applications, A James Stewart (1997), IEEE Trans on Vis & Comp Graphics

Please note that I use n here as the number of points along one dimension of a 2D raster. If you choose n = total points in the raster, algorithm complexities are O(n) and O(n(log(n))^2) respectively.

I am curious if anyone in the GRASS community has used either of them what the experience was like? Do they perform as well in practice? Is either worth considering as an upgrade to the current version of r.horizon?

[petrasovaa]

Thanks for looking into it. I don't have any experience with these, but I didn't have performance issues with the current implementation. Is this something you would be interested to implement? @metzm, any opinion on the algorithms?

[girishnand]

Thanks.
I will update as and when I get something working.