PDERL: an accurate and fast algorithm with a novel perspective on solving the old viewshed analysis problem

Viewshed analysis based on the regular grid digital elevation model (DEM) is one of the basic functions of geographic information systems. Traditional viewshed analysis algorithms are mainly carried out in a geospatial coordinate system, that create complexities when it is necessary to accurately express a perspective relationship. Moreover, it can only reduce the amount of calculation by using approximation methods, which result in decreased accuracy. We found that the error ratio of approximation algorithms is low but uneven and thereby highly likely to result in continuous area visibility errors. Despite the quality of DEM data having improved significantly in the last decade, traditional approximation algorithms have been unable to take full advantage of the improved data. Finding a fast and accurate algorithm has therefore become an urgent task. This study takes a new perspective to solve this problem by establishing a ‘proximity-direction-elevation’ (PDE) coordinate system and proposing a PDE spatial reference line (PDERL) algorithm within it. Many experiments prove that the accuracy of PDERL is the same as the R3 algorithm and its speed is about half that of the XDraw and reference plane algorithms, but much faster than R3. Due to the speed of PDERL, it is now possible and practical to do an accurate large-area DEM-based viewshed analysis.


Introduction
Digital elevation model (DEM) based viewshed analysis is one of the basic functions of many geographic information systems (GIS) (Qiang et al. 2019;Wang et al. 1996) and has broad applications. Teng et al. (1993) used viewshed analysis for topographic path planning, Winter-Livneh et al. (2012) and Bongers et al. (2012) used it for archaeological research, and Qiang et al. (2019) used it for visibility analysis of oceanic blue space. The study of viewshed analysis is based on the concept of line of sight (LOS), which refers to the line connecting an observer and a target (Franklin et al. 1994). When investigating visibility, it is necessary to know the standing positions of an observer, a target, and the intervening terrain. If any position along this line is blocked by terrain, the observer will not be able to observe the target; this phenomenon is referred to as invisibility while the converse is called visibility. Viewshed analysis refers to identifying ground areas that can be seen from a specified position within a certain range (Wang et al. 1996).
The spatial relation of viewshed analysis is very simple; however, the complexity lies in the number of calculations required. If no approximate methods are adopted, the number of calculations within traditional algorithms will be quite large. The R3 algorithm (Franklin et al. 1994) is the most direct and accurate viewshed analysis method that considers the impact of each grid line between the observer and target on the visibility. However, because the visibility of almost every target needs to be calculated individually, the computational burden is enormous, and the algorithmic complexity is O(r 3 ), which limits its practical application, especially in large-area situations. The optimization of the computational requirements has always been a research hotspot in viewshed analysis. By using approximate methods, the R2 and XDraw algorithms proposed by Franklin et al. (1994), reduced the complexity of a viewshed analysis algorithm to O(r 2 ); the reference plane algorithm proposed by Wang et al. (1996), pushed Supplementary Information The online version contains supplementary material available at https://doi.org/10.1007/s12145-020-00545-7. the speed of the single-thread viewshed analysis algorithm to its peak at the time. On the basis of these traditional foundation methods and other similar studies, some scholars went on to improve the algorithm speed by reducing the computation on a single LOS (Liu et al. 2011;Ying 2005), optimizing data exchange (Fishman et al. 2009;Haverkort et al. 2009), adopting cluster computing (Llobera et al. 2010;Song et al. 2016), or using graphics processing unit (GPU) hardware (Cauchi-Saunders and Lewis 2015; Zhao et al. 2013). Research on improving viewshed analysis algorithm speed remains an important area of study Dou et al. 2019;Li 2018).
Over the years, there have been numerous achievements in substantially increasing the speed but not the accuracy of viewshed analysis. Historically, most researchers believed that compared to the acquisition error of the DEM data itself, the calculation errors of fast approximation algorithms were relatively small, and their influence could be ignored. However, this study challenges that view. With the development of sensor technology, the resolution and precision of DEM products has been increasing, and the acquisition error of DEM data has been significantly reduced, especially in the last decade. For example, the AW3D enhanced product released by NTT RESTEC in 2014 had a resolution of 0.5 m, and the NEXTMAP ONE™ product released by Intermap in 2016 had an official elevation accuracy of 1 m ; this trend will improve with time. Traditional approximate fast algorithms cannot take advantage of the precision of modern DEM products. Furthermore, in addition to the error rate, we should also be concerned about error-point clustering.
Although the algorithm has a low error rate if the error points are clustered together it will substantially impact the results. In short, high-quality viewshed analysis not only requires higher algorithm speed, but also cannot ignore the loss of accuracy. While the XDraw and reference plane algorithms are still the two most popular algorithms for the current parallel research (e.g. Cauchi-Saunders and Lewis 2015; Dou et al. 2019;Li 2018;Song et al. 2016), their calculation errors have the greatest impact on viewshed algorithms, and most of the improvements still focus on speed while ignoring accuracy. XDraw is the least accurate algorithm among the traditional approximation methods (Kaučič and Žalik 2002), and reference plane has the same accuracy problem because of their similar perspective relationship. This study finds that although the error rate of these traditional approximation algorithms is not high; the distribution of error points is uneven. In our experiment, more than half of the cases had at least one visibility error in a contiguous area, except over plain terrain, that would seriously affect the validity of the results. Such consequences are sometimes unacceptable; for example, in military applications, where a strategy made based on an unreliable calculation could lead to unintended and deadly consequences.
Because of the difficulty of accurate perspective calculation in geographic space, and the stereotype about the low accuracy of traditional DEM, only a few studies have investigated the issue of accuracy problem of viewshed analysis. Izraelevitz (2003) improved the accuracy of XDraw by backtracking along the LOS. The accuracy of this, however, depends on the number of backtracking intersections, thus making it more like a combination of XDraw and R3. Zhi et al. (2011) proposed that the reference plane algorithm was flawed for topographic calculation in areas such as peaks and valleys, and they proposed a solution to this problem by recording the historic minimum visible elevation of each target point and comparing to the LOS height. Although its accuracy improved by 56.91%, the computation time increased by 9.34 times (Zhu et al. 2019). Zhu et al. (2019) proposed an algorithm called HiXDraw to limit the chunk distortion of XDraw by finding the most related contributing points of the target. According to their research, this method achieved a 34.96% improvement in accuracy with 4.4 times the computation time, when compared with XDraw. Their research reduced the error rate of traditional approximation algorithms, but did not eliminate it, and it was achieved at a great cost to efficiency. The studies above show that the improvement of precision is of utmost importance, especially with respect to the problem of error point aggregation. However, due to the design of the algorithm itself, it is difficult for the traditional approximation algorithms to achieve advantages in both accuracy and speed.
In addition to the development of traditional algorithms, some new solutions have emerged. Wang et al. (2015) and Fang et al. (2011) used models to replace DEMs as the basic data for viewshed analysis, and integrated rendering and calculation into a GPU. In particular, the study by Wang et al. (2015) enhanced both algorithm speed and accuracy, presenting certain advantages over traditional DEM-based algorithms. However, this does not mean that viewshed analysis based on DEM is no longer valuable. Generally, the main function of viewshed analysis is not only to mark areas visually, but also to provide basic information for other analysis algorithms. In most cases, the results do not need to be displayed, and the overhead associated with model construction, analysis, and rendering are not necessary. To reduce the pressure of computer 3D visualization for large scenes, it is typical for many systems to simplify large-scale models through levels of detail (LOD) technology, model simplification algorithms, tile block loading, and other technologies (Li et al. 2017;Wu 2019). However, model-based viewshed analysis algorithms cannot apply the above technologies and maintain the accuracy of wide-range view calculations. The process of constructing accurate models from DEM data and then rendering them in wide-range calculations is expensive, and the overhead costs of the above steps have not been reported in the relevant literature. At present, DEMs are still the main method for terrain expression. When processing largearea viewshed analysis, the DEM-based algorithm remains the simplest and most direct method.
It is critical to reduce viewshed analysis errors while maintaining the algorithm speed. In contrast to the viewshed analysis studies in recent years which mainly focus on the parallelization of traditional fast algorithms, this study focuses on finding an accurate-fast algorithm that can replace the traditional approximate-fast algorithm. The three-dimensional geographic space, in which traditional basic non-parallel algorithms are always derived (e.g. Franklin et al. 1994;Izraelevitz 2003;Qiang et al. 2019;Wu 2001;Zhi et al. 2011;Zhu et al. 2019), leads to very concise position descriptions; however, perspective problems like viewshed analysis are complex and therefore the computational load can only be reduced by approximate calculation. If we use a coordinate system suitable for perspective, we can obtain an algorithm with both speed and accuracy. According to this principle, a coordinate system of 'proximity-direction-elevation' (PDE coordinate system) is established in this study. Its advantage is that the spatial linear relationship of lines in geographic space can be transformed into the linear relationship of direction, elevation, and proximity, which lays a foundation for fast and accurate viewshed analysis. Furthermore, a PDE spatial reference line algorithm (PDERL) is proposed based on the above coordinate system. After numerous experiments, we can prove that the speed of this algorithm is close to that of XDraw and reference plane algorithms and can therefore meet the speed requirements of practical applications. Moreover, this algorithm is precise, and its results are the same as the R3 algorithm, thus solving the problem of error point aggregation in the traditional approximate algorithms. The algorithm proposed in this study can alleviate the long-standing trade-off between speed and accuracy of DEM-based viewshed analysis algorithms.

Materials
The experimental materials in this study consists of three DEMs from ASTER GDEM with 30 m grid spacing, in the TIFF format, whose serial numbers are N34E114, N41E119, and N28E097. The locations and topographic features of the three DEMs are shown in Table 1, and the height images of the DEM data are shown in Fig. 1. The three DEMs have a total area of approximately 36,000 km 2 , and are broadly representative-covering typical landforms, including plains, hilly areas, and mountains-to ensure that the analysis results are universal. The algorithm in this study and other comparison algorithms are used to compare viewshed analyses based on these materials.

Fundamental basis
The visibility of the target is determined by the perspective relationship between the target and an obstruction. If there is an occlusion on the terrain, the following three conditions must be met simultaneously.
(1) The obstruction is in front of the target.
(2) The obstruction and the target are in the same direction from the observer. In this study, the PDE coordinate system is proposed to directly quantify the three conditions above. The description of the object in this space is transformed from the traditional geometric position coordinates into three conceptual dimensional coordinates relative to the observer: proximity (p), direction (d), and elevation (e). The PDE coordinate system is a right-handed coordinate system with p,d, and e as the X-axis, Y-axis, and Z-axis, respectively. If the two targets are in the same direction (d is the same) and the larger the value of e, the greater the elevation angle of the target LOS. Moreover, and the higher the value of p, the closer the target to the observer.
A DEM grid point is a spatial point with definite coordinates and elevation. Assuming that the grid line between the two points adjacent to the coordinates on the DEM can express the corresponding terrain, the continuous polyline formed by the DEM points of a row/column can express the terrain of the corresponding geographic space of the row/column, which is referred to as the latitude/longitude terrain line in this study. By definition of the spatial characteristics of PDE, any point located between and colinear with two points in geographic space is still located between and colinear with them in PDE space, so the polyline formed after terrain line is converted to PDE space can also represent any point on the original terrain line. When carefully dividing the calculation regions, any terrain line converted to its corresponding PDE space will be in a plane with a constant proximity value and will also be a continuous polyline. The viewshed problem on a DEM can be abstractly discussed as the occlusion of the terrain line close to the observer (with larger proximity) to the terrain line far away from the observer (with smaller proximity). A latitude/longitude terrain line has an independent effect on the occlusion of grid points, but a grid point can be seen by the observer only if it is not shielded by both. Therefore, based on the position of the observer, this study divides the latitude and longitude terrain lines in the calculation region into four regions, namely Regions I, II, III, and IV, as shown in Fig. 2. Then each region is converted into the corresponding PDE space for calculation, and the visibility of any grid point will be judged by the latitude (Regions I and II) and longitude (Regions III and IV) terrain lines. The result can be obtained by combining the two. As this algorithm is based on the PDE coordinate system, this study introduces the definition and properties of the PDE coordinate system first, and then provides the derivation of PDERL in detail.

Definition of PDE coordinate system
The PDE coordinate system is defined based on a topocentric coordinate system. We define the topocentric coordinate system for each of the four regions. The topocentric coordinate system for Region I is defined as follows: take the point of view of the observer as the origin, the north direction of the DEM vertical grid line as the Y-axis, the east direction of the DEM horizontal grid line as the X-axis, and the positive direction of a normal vector on the reference ellipsoid of the observer as the Z-axis to establish the right-hand spatial rectangular coordinate system. While taking the Z-axis of the topocentric coordinate system in Region I as the axis of rotation, we can rotate the coordinate system clockwise to 180°, 270°, and 90°, to obtain other topocentric coordinate systems in Region II, Region III, and Region IV, respectively. The geographic location of grid nodes in DEM files is generally recorded by the coordinates of longitude and latitude or a specific kind of map projection coordinates (such as Universal Transverse Mercator (UTM) projection), which record positions in space independent of the observer. When the coordinates of grid nodes are converted into the above topocentric coordinate system, they will be in relation to the relative spatial position of the observer. This study defines three concepts-proximity, direction, and elevation-based on the topocentric coordinate system and uses them to quantify the three conditions of occlusion. The space constituted by the quantized value of these three concepts is defined as PDE space. This study introduces the definition of PDE space using point A as an example in Fig. 3a, and the same point transformed to PDE space as shown in Fig. 3b. We defined four topocentric coordinate systems that can be overlapped by the rotation described above. That said, the steps to convert coordinates from any topocentric coordinate system to PDE space coordinate system are the same. The definition below in this section will refer to the topocentric coordinate systems in Region I as an example. It should be emphasized that all coordinates (x, y and z) in the definitions of proximity, direction, and elevation below are the values in the topocentric coordinate system, not the coordinates of latitude and longitude or the map projection coordinates recorded by DEM itself.
The proximity p is defined as the reciprocal of the X-axis value of the investigation point in the topocentric coordinate system, and the definition is expressed as formula (1). The domain of p is(−∞, 0) ∪ (0, +∞). When two objects are in the same direction, the object with a larger proximity value (i.e. which is closer to the observer) may block the object with a smaller proximity value (i.e. which is farther away from the observer). This quantifies the occlusion condition (1).
The direction d is defined as the tangent value of direction angle α in the topocentric coordinate system of the investigated point, and is defined by formula (2). The domain of d is(−∞, +∞). When the d values of the two points are different, they are in different directions from the observer's perspective. This quantifies the occlusion condition (2).
The elevation e is the product of the relative elevation z and proximity p of the observation point and is defined by formula  (3). The domain of e is(−∞, +∞). As in Fig. 3a, the meaning of e in geographic space is the tangent value of elevation γ, after point A is shifted along the Y-axis toA ′ . The slope (s) is the tangent value of elevation φ in the topocentric coordinate system; the higher the value of s, the higher the elevation angle of the LOS. Meanwhile, e and s conform to the relation formula (4). Therefore, when d is the same, the larger the value of e, the larger will be the value of s in the LOS. This quantifies the occlusion condition (3).
The PDE coordinate system is defined as follows: take p as the X-axis, d as the Y-axis, and e as the Z-axis in a right-hand coordinate system. Combined with some mathematical principles, it is easy to prove that for any straight line in the topocentric coordinate system that is not in the Z-O-Y plane, when converted to the PDE space, the p, d, and e of the line are linearly related to each other, indicating that it remains a straight line in PDE space. This property lays a theoretical foundation for fast accurate calculation of the viewshed algorithm in this study.

The essence of viewshed analysis in PDE space
The transformation from the topocentric coordinate system to PDE space coordinate system can be visualized as a complex perspective projection with an observer transformed into an orthographic projection where the observer is set at an infinity distant point on the P-axis and is looking in the opposite direction of the P-axis. When a terrain line has been converted to PDE space, it will be a perspective drawing of the terrain, and we call it the investigation line. When calculating the occlusion of an investigation line with a larger p to an investigation line with a smaller p, the former may block the later. The investigation line with a larger p is a reference when calculating the occlusion of any line with a smaller p, so we can call the former a reference line to the latter. Like the characteristic of orthogonal projection, we can judge visibility by the method of "overlay" of the investigation lines.
Taking Region I as an example, we convert the first terrain line C 1 to PDE space to form the initial reference lineL R , as shown in Fig. 4, and convert the second terrain line C 2 to PDE space to constitute investigation lineL 2 . Next, we shift L R along the proximity axis to the plane where L 2 is, and divide the terrain line into several sections based on the intersection points between L R andL 2 . In fact, at the intersection of the two terrain lines, the visibility state of C 2 changes. The sections of the investigation line with a lower e will be covered and the grid points located on these sections will also be covered. Finally, we take the sections with higher e of L R and L 2 to form a new terrain line L 0 R , which in turn will be used as the reference line when calculating the terrain line C 3 in the next steps. Following these steps, we started from the first two terrain lines closest to the observer in one region and calculated the terrain line one by one according to the above steps until it was completed. The calculation steps of one region are shown in Fig. 5.

Expression of terrain line in PDE space
A terrain line is a spatial polyline constructed by a number of continuous DEM grid points. According to PDE spatial properties, for any segment of a terrain line there must be appropriate constants a and b that make the line segment conform to a b formula (5) in PDE space (see Supplementary Material 1 for the derivation). As p = 1/x on each terrain line is constant, the terrain line must be a continuous polyline located in a PDE plane with constant proximity after being converted into PDE space. In this study, the line segment converted to PDE space is expressed as L(a, b| d l , d r ), where d l and d r are the left and right boundaries of the defined domain in the direction of the line segment and a, b are the parameters of L. If the topocentric coordinates of two adjacent grid points are known to be A(x, y A , z A ) and B(x, y B , z B ), and u is marked as the reciprocal of the grid width of two adjacent grid points, then the value of constant a can be expressed as formula (6) (see Supplementary Material 2 for the derivation).
In fact, b is not necessary to describe a terrain line. By recording e of the first segment, a of each segment, and the ending d of each segment, we can express and record the whole terrain line in PDE space.

Discriminant of visibility
As shown in Fig. 6, the difference ΔE C by e C , e F is defined as formula (7), corresponding to the vector in Fig. 6. If r C is used to mark the visibility of point C, the value of 0 indicates invisibility, and the value of 1 indicates visibility; then, we can easily obtain the discriminant of visibility as formula (8).

Construction and update of reference line
The calculation of the visible segment on the investigation line and the update of the reference line are a cumulative process of ΔE. As shown in Fig. 7, the current investigation line is labelled L C and the current reference line is labelled L R which has been moved to the plane where L C is located. The  direction values of all the inflection points of L R and L C are taken as the segmentation points with which the two terrain lines are divided into several segments. The direction at the end of segment i is labelled d i , and the values of parameter a on L R and L C in this segment are a C i and a R i . The value of ΔE at the end of this segment is ΔE i , while at the end of the previous segment it is ΔE i -1 . Due to the linear characteristics, it is easy to obtain formula (9), the cumulative formula for ΔE. If the end of a section is a grid point on the investigation line, the visibility of the grid point can be obtained by a combination of formulas (8) and (9). If ΔE i − 1 and ΔE i have opposite positive and negative values, then L R and L C have an intersection point in segment i, and the direction of the intersection point (d ′ ) can be solved by formula (10). By using formulas (7)-(10), we can follow the steps in Fig. 8 to calculate the visibility of one investigation line and update the current reference line.
When we follow the last step of Fig. 8, the identified intersection points (generally not grid points) were added to the new reference line as inflection points when we selected the sections with higher e to update the reference line. Thanks to the precise representation of the reference lines by these nongrid points, the accuracy of the PDERL algorithm results will be comparable to those of the R3 algorithm. The process of updating the reference line is the main computational process of PDERL, and the formulas (8)-(10) will be the most utilized in this process. We derived these formulas in the PDE coordinate system; however, it is not necessary for an independent step to convert DEM from its recorded coordinate system into the PDE coordinate system in the actual calculation. As all the above formulas can be converted into equivalent mathematical formulas expressed by geographic space coordinates through the definitions, when we use them for viewshed analysis, it already contains the coordinate transformation process.

Consolidation viewshed results of the four regions
Grid points can be judged to be visible only when they are determined to be visible by the longitude terrain lines (in Region I/II) and the latitude terrain lines (in Region III/IV) simultaneously. In this study, a 0-1 matrix is used to record Fig. 7 Relationship between ΔE i and ΔE i -1 . The figure shows a PDE plane where L C is located, and L R has been translated from a plane with a bigger proximity to this plane. L R is the reference line to L C . The values of parameter a of AB and MN are a C i and a R i . It is evident from the figure that Þ , which is the theoretical basis of the cumulative formula for ΔE Current reference line LR , and the investigation line LC .
Tak in g the direction valu es of all inflection points of LR and LC as segmentation points, LR and LC can be divided into n pairs of line segments t hat contain both of them. The segmentation point i on LR and LC is marked as Pi R and Pi C respectively.
Use formula (7) to calculate the ΔE of P0 R and P0 C ; use formula (8) to judge the visibility state of P0 C , and mark the result as Vlast; record Vlast as the visibility state of the first grid point on LC ; set i=1; Use formula (9) to calculate the ΔE of Pi R and Pi C ; use formula (8) to judge the visibility state of Pi C , and mark the result as Vcurrent . End the calculation of LC and the updating process of LR .
If the Pi C is a grid point, record the visual state. Fig. 8 Steps of visibility calculation of one general investigation line viewshed analysis results, which have the same dimensions as the viewshed DEM, and its default value is 0. The results of viewshed analysis for each grid point are stored in the corresponding element with the same row and column number. The value of 0 represents invisibility, and 1 represents visibility. The results of Regions I, II, III, and IV are recorded in the visibility matrix V 1 , V 2 , V 3 , and V 4 . The final visibility matrix is obtained using formula (11).

Experiments and results
To verify the effectiveness of the algorithm, we designed four contrast experiments on PDERL, R3, XDraw, and reference plane algorithms to ascertain accuracy, aggregation degree of error points, computational complexity, and speed. As a classical viewshed analysis algorithm, R3 considers all grid lines in LOS as they affect the visibility. The results of R3 are taken as the standard because of their high accuracy and reliability and are used to evaluate the error points of the results of the other algorithms. This study implements all the compared algorithms with C# programming in a back-end web service, while developing a front-end interactive platform with Cesium. These algorithms run on a Mi computer 181,501-AB with a Windows 10 operating system.

Experiment 1: Accuracy
This experiment began at an altitude of 1 m above the ground. For every 25 m increase in height, the observation location and analysis radius were randomly selected (radius: 270-10,680 m). The experiments were conducted 50 times for each of three DEM blocks. For each experiment, PDERL, XDraw, reference plane, and R3 were used to calculate the same region, and the results of each algorithm were compared with the results of R3 to analyze the corresponding errors. The statistical results are shown as polynomial fitting curves of the experimental data in Fig. 9. By comparing the upper limit of error rate, average error rate, and median error rate of each algorithm, it can be observed that: (1) The results of PDERL in all experiments are consistent with R3; (2) The error of the XDraw/reference plane algorithm is the largest near the ground and decreases gradually with the increase in observation height; (3) The XDraw/reference plane algorithms have a higher error rate near the ground in the flat terrain area than in the undulating terrain area, but it declines rapidly with the height from the ground.

Experiment 2: Aggregation of error points
The PDERL algorithm has no error points; accordingly, there is no possibility of aggregating error points as in most experimental results of XDraw and reference plane algorithms. For example, Fig. 10 shows the viewshed analysis results calculated by PDERL and XDraw algorithms in an experiment. The coordinate of the observation center is 119.8219°E, 41.7602°N , 5 m above the ground, and the observation radius is 1 km. As shown in the figure, there are some error point aggregations in the XDraw result while there are no error points in the PDERL result. In fact, the error rate of XDraw algorithm is not high in this experiment, only 3.1221%, but such a low error rate still cannot prevent the problem of many error points becoming aggregated.
To quantify the degree of aggregation of the error points in the results, this study defines the concept of an 'x-neighborhood area error (x-NE)' and an x-NE logical flag (x-NEF) corresponding to whether one result consists of at least one x-NE (1) or not (0). An x-NE means the area where 90% of the center grid point and the grid points in its x-grid distances around (as shown in Fig. 11) have the same kind of visible error. For the x-NEF, first select an algorithm to perform a viewshed analysis and then check whether the analysis result consists of at least one x-NE. The value of the x-NEF is 1 for  Fig. 9 Error rate statistics of the algorithms on three terrains. The left/middle/right subfigure is the error rate of limit/average/mean square of PDERL, XDraw and reference plane algorithms sorted by height above the ground. Since the XDraw and reference plane algorithms have the same error condition, they are drawn on the same line yes or 0 for no, instead of a continuous value, and its main function is to count how many results are affected by x-NEs by summing the x-NEFs. The DEM resolution tested in this study is 30 m, and the size of the area or grid points involved in various types of x-NEs are shown in Table 2.
In this experiment, heights from 2 to 52 m above the ground (at 2 m intervals) for three DEM sections (covering 120.5604 km 2 each) with a fixed radius (5.49 km) and randomly selected location were computed 500 times with R3, PDERL, and XDraw. This comparative experiment was carried out 12,500 times on each of the three DEMs and 1500 times for each height. Considering the large number of experiments undertaken, the results were representative to some extent. In this study, there was no error point in the results of PDERL, while there were many error points in the results of Xdraw. The x-NEFs of the results of Xdraw algorithm were counted as shown in Fig. 12 and Table 3 according to their height order and DEM file. As can be seen from the relevant figures and tables, (a) in mountainous (N28E097) and hilly (N41E119) areas, more than half of the calculation results have 1-NEs; (b) the larger the range of x, the lower the probability of x-NEs. Furthermore, in the calculation results of the mountainous area (N28E097) the occurrence frequency of 10-NE is still 0.15%, which means that in this experiment 0.15% of N28E097 calculation results experienced at least 396 consecutive grid-point errors, involving a continuous error area of up to 396,900 m 2 .

Experiment 3: Computational complexity
The PDERL and R3 algorithms are both exact algorithms; however, the complexity of R3 is O(r 3 ) and the theoretical complexity of PDERL is the same as XDraw and reference plane, which is O(r 2 ). Theoretically, PDERL should   be much faster than R3 and closer to XDraw/reference plane. To compare their computational complexities, an experiment was conducted on the DEM N41E119 at a height of 2 m from the ground, radius from 300 m to 4800 m (at 30 m intervals), computed 11 times with R3, PDERL, XDraw, and reference plane on randomly selected locations. In this study, the computation times of each algorithm was plotted statistically in order of radius size (Fig. 13). The experimental results confirm our assumption on the degree of computational complexity of the considered algorithms.

Experiment 4: Speed
To compare the speed of PDERL to that of XDraw, reference plane, and other recent improvements in accuracy of The sum of x-NEFs on hilly terrain (N41E119); (c) The sum of x-NEFs on mountainous terrain (N28E097). Site height is the height above the ground Note that the distribution ratio is obtained by dividing the sum of x-NEFs on different DEM by the total number of experiments on the corresponding DEM, respectively. It is an indicator used to measure the probability that a kind of x-NE will affect the result of a viewshed analysis.
algorithms, for example as in Zhu et al. (2019), this study conducted 6 experiments with each of the three algorithms, on the three DEMs, at heights from 1 m to 5000 m from the ground (in 1 m intervals), with a radius of 10.77 km, at randomly selected locations, totaling 15,000 comparison experiments. We used a concept called time consumption ratio (R), similar to Zhu's approach (Zhu et al. 2019), to evaluate the speed difference between PDERL, XDraw, and reference plane. Time consumption of PDERL, XDraw, and reference plane were defined as T PDERL , T XDraw , and T rp , respectively. They were compared in an experiment through their time consumption ratio between XDraw (R XDraw ) and reference plane (R rp ) in n experiments with the same height as defined by formula (12).
After fitting with a 6th degree polynomial, the results of R XDraw and R rp in the order of height in different terrains are shown in Fig. 14 below. The speed performance of the PDERL algorithm is not very stable, particularly at a low height, but the overall difference in speed is not large, a little faster than half the speed of XDraw, and a little slower than half the speed of reference plane. Although we did not implement the method in Zhu et al. (2019) we made a comparison based on their research. Their algorithm is 4.4 times as time-consuming as XDraw with limited accuracy improvement (Zhu et al. 2019), while PDERL is at most 2.25 times as time-consuming as XDraw with completely accurate results.   Fig. 14 The time consumption ratio of PDERL to XDraw and reference plane. The left/right subfigure shows the time consumption ratio of PDERL to XDraw and reference plane at different heights from the ground. While the time consumption of this experiment is relatively stable at high altitudes, it is approximately 185% of XDraw and 215% of the reference plane, and the maximum value of the experiment is not more than 225% of XDraw and 255% of reference plane

Compliance with ethical standards
Conflicts of interest/competing interests No potential conflict of interest was reported by the authors.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.