Algorithms for viewshed calculations have been reviewed on and off over the past several weeks. Based on this background reading, I have been developing the basic concept of an algorithm that only refers to the neighbour points in calculating viewability from a location. I wanted to test out the concept before looking at extending it to 3D.
The general concept was taken from the below reference
- A Fast Algorithm for Approximate Viewshed Computation ( https://pdfs.semanticscholar.org/63c1/857866341dbf79158fdbfa4c9f4046a11e3d.pdf )
Develop base algorithm for viewshed calculation in order to inform planning for 3D implementation
- Test scene created in Unity using cubes in a side on view (2D)
- Extracted data from scene (cube heights in y axis) to provide initial conditions for algorithm
- Implementation of viewshed algorithm
- Added observer height and observation tolerance (max viewshed height to be considered out of Line of Sight – eg. If a point is 1cm below view it can still be considered ‘in view’)
- Implemented visualisation in Unity and randomised location demo
The algorithm in a linear sense is basically as follows:
- Given a set of heights based on an x position (eg. x = 0, y = 1.3, x = 1, y = 1.6 etc)
- Create an array of heights with the index of each array element referring to the x position
- Create an equal sized array of Viewshed heights (Corresponding to the height something must be at to be visible from the viewing location
- Given a viewer location (x position)
- Set the viewshed height of that index and those either side of it to 0 (ie. you can always see the point you are on and the neighbour points
- Moving to the right (increasing x position value)
- Calculate gradient at previous position = rise / run
- RISE = (height of previous position + viewshed height at previous position) -height at viewers location
- RUN = Absolute value of (viewerLocation index – previous location index)
- Given that gradient
- Minimum height to be seen at current index = height at previous index + viewshed height at previous index + gradient
- (The above is due to us moving ONE unit at a time on the X axis)
- IF height at current index is greater than minimum height to be seen at the current index, set viewshed height at current index to 0 (We can be directly seen)
- OTHERWISE set the viewshed height at current index to minimum height to be seen at current index – height at current index
- This works well as can be seen in the gif below
In 3D the viewer point will be in an arbitrary position in an X Y plane. My initial thoughts are the exact same (linear) algorithm can be used in the cardinal directions (ie. ‘North/South’ being directly along the Y axis, ‘East/West’ being directly along the X axis from the viewer position only). A very slight modification can also be used for the NW/NE/SW/SE directions as well.
The algorithm will need to be amended for all other cells however; It will not be possible to directly check a previous location’s height and viewshed height. This is not a significant complication though; it is a relatively straightforward process to interpolate between the two points closest to where the ‘previous neighbour’ would be (along the Y axis for locations in the North and South quarters and along the X axis for locations in the East and West quarters). This will involve a weighted average between the two points with the weighting based on the relative distances on where the exact neighbour should be (eg. If the neighbour location would be 4.2 on the x axis, the values from the point at 4 would be weighted at 80% compared to the values at the point at 5 on the x axis with a weighting of 20%)
One significant outcome of the above is, once the cells in black above are calculated with the simple algorithm, each of the 8 subsections created can have their values calculated independently of the other sections. This opens up potential for a multi threaded approach which should be relatively straight forward (each segment is separated from the others and has no interaction with other elements) and will be an increase in efficiency.
Original: 4 hours (24 hours running total)
From update below: 5 hours (29 hours running total)
Work continued on the above throughout the day. The basic linear viewshed class was refined and polished to a standard I am very happy with. The general algorithm for 3D was fleshed out and it was implemented in a Viewshed3D class. This is untested as yet as I will need to alter the visualisation scripts somewhat to effectively see what is happening ‘under the hood’. I ended up just breaking the space into quadrants (separated by the X and Y lines from the viewer location) and this provides a much easier implementation.
Work was also started on drafting the final report; I will start fleshing parts out now (especially the background/planning elements).