Efficient plotting the functions with discontinuities based on combined sampling Efficient plotting the functions with discontinuities based on combined sampling Tomáš Bayer Department of Applied Geoinformatics and Cartography, Faculty of Science, Charles University, Czech Republic bayertom@natur.cuni.cz Abstract. This article presents a new algorithm for interval plotting of the function y = f(x) based on combined sampling. The proposed method synthesizes the uniform and adaptive sampling approaches and provides a more compact and efficient function representation. Dur- ing the combined sampling, the polygonal approximation with a given threshold α between the adjacent segments is constructed. The automated detection and treatment of the disconti- nuities based on the LR criterion are involved. Two implementations, the recursive-based and stack-based, are introduced. Finally, several tests of the proposed algorithms for the different functions involving the discontinuities and several map projection graticules are presented. The proposed method may be used for more efficient sampling the curves (map projection graticules, contour lines, or buffers) in geoinformatics. Keywords: function; adaptive sampling; combined sampling; recursive approach; stack; dis- continuity; polygonal approximation; visualization; map projection; plotting; GIS; Octave; Mathematica. 1. Introduction A function y = f(x) on interval Ω = [a,b] may have different form. For efficient plotting, its polygonal approximation needs to be constructed. A current approach concentrated on uniform sampling with the step δx may not be sufficient. Despite its popularity, the equally spaced points cannot describe its course without errors; the problems of undersampling or oversampling are common. Adaptive sampling brings several benefits, it adapts to a different curvature of the function, reduces the amount of redundant data and provides a natural and smooth plot of the function without the jumps and breaks. This technique is popular in computer graphics; recall the well-known deCasteljau or Chaikin’s algorithms for the curve approximation. Combining the uniform and adaptive sampling approaches their advantages can be synthesized. The resulted representation is more compact, efficient and smooth. A function may contain points of discontinuities that need to be detected. A subdivision of the given interval Ω, to the set of disjoint subintervals Ωgk without the internal singular- ities, containing only “good” data needs to be undertaken. In other words, the polygonal approximation of the function needs to be split into the several disjoint parts. The pro- posed method works by the requirements mentioned above. A broad set of the singularities can be detected and treated using the multiple criteria. Subsequently, for each interval Ωgk, a polygonal approximation of f(x) is constructed using combined sampling. Since there are many sophisticated solutions built-in the high-end systems (Mathematica, Maple), our simple algorithm based on the recursive approach is efficient and easy-to-implement. This paper is organized as follows. In Section 3, the combined sampling technique for 1D functions is presented. Section 4 deals with the detection of singularities, Section 5 describes the combined sampling of the discontinuous functions. In Section 6, the combined sampling Geoinformatics FCE CTU 17(2), 2018, doi:10.14311/gi.17.2.2 9 http://orcid.org/0000-0001-6279-1926 http://dx.doi.org/10.14311/gi.17.2.2 http://creativecommons.org/licenses/by/4.0/ T. Bayer: Efficient plotting the functions with discontinuities technique is tested on the set of several functions. Subsequently, its behavior on four map projections is analyzed. 2. Related Work There are several strategies for plotting the function y = f(x) on interval Ω = [a,b]. The naive approach based on sampling of f in a fixed amount of the equally spaced points is described in [20]. The simple functions suffer from oversampling, while the oscillating curves are under-sampled; these issues are mentioned in [14]. Another approach based on the interval constraint plot constructing a hull of the curve was described in [6], [13], [20]. The automated detection of a useful domain and a range of the function is mentioned in [41]; the generalized interval arithmetic approach is described in [40]. A significant refinement is represented by adaptive sampling providing a higher sampling density in the higher-curvature regions. The are several algorithms for the curve interpo- lation preserving the speed, for example: [37], [42], [43]. The adaptive feed rate technique is described in [44]. An early implementation in the Mathematica software is presented in [39]. By reducing data, these methods are very efficient for the curve plotting. The polygo- nal approximation of the parametric curve based on adaptive sampling is mentioned in the several papers. The refinement criteria, as well as the recursive approach, are discussed in [15]. An approximation by the polygonal curves is described in [7], the robust method for the geometric and spatial approximation of the implicit curves can be found in [27], [10], the affine arithmetic working in the triangulated models in [32]. However, the map projections are never defined by the implicit equations. Similar approaches can be used for graph drawing [21]. Other techniques based on the approximation by the breakpoints can be found in many papers: [33], [9], [3]; these approaches are used for the polygonal approximation of the closed curves and applied in computer vision. 3. Combined sampling In this section, the proposed combined sampling technique providing the polygonal approxi- mation of the parametric curve involving the discontinuities will be presented. The modified method will be used for the function f(x) reconstruction and plot. Based on the ideas of splitting the domain into the subintervals without the discontinuities, it represents a typical problem solvable by the recursive approach. 3.1. Polygonal approximation of the curve Let y = f(x), M ⊂ R, f : M → R be a function of a real variable x, and the set M represents the domain of the function f. Let Ω = [a,b], a ∈ M, b ∈ M be the subdomain inside which the polynomial approximation pi = (xi,f(xi)), 1 ≤ i ≤ n of the curve f is constructed, where x1 = a < x1 < ... < xn = b. This approach leads to a discrete reconstruction of f from the set of sampled points. The behavior of f should be reconstructed concerning its curvature. The classical approach based on uniform sampling from the equidistant points xi with the step δ, where δ = xi+1−xi, provides a good approximation only if δ → 0. For the straight parts of the curve, many Geoinformatics FCE CTU 17(2), 2018 10 T. Bayer: Efficient plotting the functions with discontinuities -4 -3 -2 -1 0 x -2 -1 0 1 2 f( x ) AS: α =1°, 83 sampled points -4 -3 -2 -1 0 x -2 -1 0 1 2 f( x ) AS: α =5°, 23 sampled points -4 -3 -2 -1 0 x -2 -1 0 1 2 f( x ) AS: α =10°, 11 sampled points -4 -3 -2 -1 0 x -2 -1 0 1 2 f( x ) US: α =1°, 629 sampled points -4 -3 -2 -1 0 x -2 -1 0 1 2 f( x ) US: α =5°, 117 sampled points -4 -3 -2 -1 0 x -2 -1 0 1 2 f( x ) US: α =10°, 63 sampled points Figure 3.1: Adaptive (AS) and uniform (US) sampling of the meridian of the longitude λ = −180◦ in the Sanson projection for angles αi = 1◦, 5◦, 10◦. almost colinear segments are constructed; too much redundant data is generated. Conversely, for larger δ, the shape of the function in the high-curvature areas may not be captured adequately, which is discussed in [15]. In general, the main disadvantages of uniform sampling, the problems of undersampling or oversampling, are referred. For a current density of the sample, the equally spaced points cannot describe the function course without errors. 3.2. Combined sampling technique By avoiding colinear segments as well as a better adaptation to the different curvature, adap- tive sampling respects the behavior of the function more naturally. It leads to the more compact data representation of the curve while its shape, as well as its aesthetic look, are preserved. A comparison of adaptive and uniform sampling for the meridian curve is illus- trated in Fig. 3.1. Uniform sampling requires more data to maintain the same curvature represented by the angle αi between the adjacent segments (pi−1,pi) and (pi,pi+1). The dif- ference increases depending on the curvature. In general, for adaptive sampling, the required amount of points is about one order less. Unfortunately, two issues are referred in [15]: 1. For specific functions, some narrow subintervals in the early iterations may be skipped and stay unprocessed. 2. For some periodic functions, a refinement based on the iterative subdivision into the segments of the same length may not be successful, if the function values at these points are equal. To avoid the problems in the first case, it is natural to take advantage of both methods Geoinformatics FCE CTU 17(2), 2018 11 T. Bayer: Efficient plotting the functions with discontinuities -2 -1 0 1 2 x 0 0.5 1 1.5 2 2.5 3 3.5 4 f( x ) AS: d=1 -2 -1 0 1 2 x 0 0.5 1 1.5 2 2.5 3 3.5 4 f( x ) AS: d=2 -2 -1 0 1 2 x 0 0.5 1 1.5 2 2.5 3 3.5 4 f( x ) AS: d=3 Figure 3.2: Illustration of the proposed algorithm: an adaptive sampling of the function y = x2, x ∈ [−2, 2] for the depths of recursion d = 1, 2, 3. and propose the combined sampling method. Uniform sampling represents the initial step, later steps refining the curve approximation are provided by adaptive sampling. The second problem may overcome adding the partial randomness to the generated segments. Refinement criteria. An important role is played by the refinement criteria smoothing the polyline [15]. Suppose (pi−1,pi,pi+1) to be three consecutive sampled points of the curve. The primary criterion is represented by the angular difference of both segments αi = α(pi−1,pi,pi+2) = arccos |u ·v| ‖u‖‖v‖ ±π, u = pi+1 −pi, v = pi−1 −pi. Different criteria can provide the analogous results: recall the distance of pi from pi−1,pi+1, or, the local length ratio [29]. Combined sampling algorithm without the singularities The combined sampling algorithm is based on the idea of the hierarchical reconstruction of the curve shape, which follows the recursive approach with the multiple calls mixing the uniform and adaptive techniques. Unlike a simple algorithm discussed in [15], the proposed method can handle the singularities and discontinuities of f and requires a lower recursion depth. Initially, sampling without the treatment of singularities will be presented. Suppose d to be the current depth of the recursion, d to be the minimum and d to be the maximum recursion depth. Combined sampling returns a polynomial approximation of the curve by the refinement criteria α and the recursion depth d. Our algorithm combines the uniform and adaptive sampling techniques and subdivides Ω into a specified number of the disjoint subintervals Ωk of the similar size during each recursive step, where k = 4. Hence, Ω is split into the approximate quarters with the randomly shifting borders, which are not a direct multiple of 0.25. Initially, if d ≤ d the interval is subdivided into four subintervals regardless of α; four new segments of the polygonal approximation are created. Subsequently, when d > d, between Geoinformatics FCE CTU 17(2), 2018 12 T. Bayer: Efficient plotting the functions with discontinuities Algorithm 1 Combined sampling, the initial phase. 1: function csInit(f,L,a,b,d,d,d,�,α) 2: L = ∅ 3: ya = f(a),yb = f(b) 4: if discontinuity in a then 5: throw SingularityException (a) 6: if discontinuity in b then 7: throw SingularityException (b) 8: L ← Point(a,ya) 9: as(f,L,a,b,ya,yb,d,d,d,�,α) 10: L ← Point(b,yb) each pair of four new consecutive segments, the refinement criterion αk, k = 1, ..., 3, is eval- uated and compared to α. If αk > a and b−a > ε, two adjacent segments are created. The interval is subdivided into 2-4 new until the visually “smooth” polynomial approximation of the curve is obtained. In other words, uniform sampling is followed by adaptive sampling refining the properties of the insufficiently estimated segments. Let L = {pi}ni=1 be the polynomial approximation of f(x) and Ω = [a,b] a subdomain. The algorithm may be summarized as follows: 1. The initial phase Let L = ∅ be the empty set. Compute ya = f(a) and yb = f(b). If a singularity in a or b is detected, throw the exception. Add the initial vertex to L: L ← pa, where pa = (xa,ya). Set the recursion depth d = 1. 2. The recursive step Enter the recursive procedure and do the following substeps: (a) If d > d or b−a < ε stop the recursive procedure and go to Step 3. (b) For a given Ω = [a,b], the interval is split by the three points x1 = a + 1 2 r1(b−a), x2 = a + r2(b−a), x3 = a + 3 2 r3(b−a), into the approximate quarters Ω1 = [a,x1], Ω2 = [x1,x2] Ω3 = [x2,x3], Ω4 = [x3,b], where r1, r2, r3 are the random numbers inside the interval [0.45, 0.55]. This step prevents a situation, when f(a) = f(b) = f(x1) = f(x2) = f(x3), but their inter- mediate points do not held this condition; it is typical for some periodic functions (for example, if y = sin 2x, and Ω = [0, 2π]). (c) If a singularity in x1, x2, or x3 occurs, throw the new exception with the argument indicating the singularity. (d) Evaluate the function values y1 = f(x1), y2 = f(x2), and y3 = f(x3) at new vertices p1, p2, p3. Geoinformatics FCE CTU 17(2), 2018 13 T. Bayer: Efficient plotting the functions with discontinuities Algorithm 2 Combined sampling, the recursive phase. 1: function cs(f,L,a,b,ya,yb,d,d,d,α) 2: if d > d∨ (b−a < ε) then 3: return 4: r1 = rand(0.45, 0.55), r2 = rand(0.45, 0.55), r3 = rand(0.45, 0.55) 5: x1 = a + 12r1(b−a), x2 = a + r2(b−a), x3 = a + 3 2r3(b−a) 6: if discontinuity in xi then 7: throw SingularityException (xi), i = 1, 2, 3 8: y1 = f(x1), y2 = f(x2), y3 = f(x3) 9: pa = Point(a,ya), pb = Point(b,yb), pi = Point(xi,yi), i = 1, 2, 3 10: α1 = α(pa,p1,p2), α2 = α(p1,p2,p3), α3 = α(p2,p3,pb) 11: if (α1 > α) ∨ (d <= d) then 12: cs(f,L,a,x1,ya,y1,d + 1,d,d,α); 13: L ← Point(x1,y1) 14: if (α1 > α) ∨ (α2 > α) ∨ (d <= d) then 15: cs(f,L,x1,x2,y1,y2,d + 1,d,d,α); 16: L ← Point(x2,y2) 17: if (α2 > α) ∨ (α3 > α) ∨ (d <= d) then 18: cs(f,L,x2,x3,y3,y3,d + 1,d,d,α); 19: L ← Point(x3,y3) 20: if (α3 > α) ∨ (d <= d) then 21: cs(f,L,x3,xb,y3,yb,d + 1,d,d,α); (e) Check the refinement criteria α1 = α(pa,p1,p2), α2 = α(p1,p2,p3), α3 = α(p2,p3,pb), and the recursive depth d. When the curve is not sufficiently smooth, or d ≤ d, it needs to be refined. For d ≤ d, this step begins with uniform sampling; for d > d it transforms to adaptive sampling. (f) If α1 > α, call the recursive procedure with the increased depth d = d + 1 for the interval [a,x1). (g) Add new point p1 to the polynomial approximation of f(x): L ← p1. (h) If α1 > α∨α2 > α, call the recursive procedure with the increased depth d = d+ 1 for the interval (x1,x2). (i) Add new point p2 to the polynomial approximation of f(x): L ← p2. (j) If α2 > α∨α3 > α, call the recursive procedure with the increased depth d = d+ 1 for the interval (x2,x3). (k) Add new point p3 to the polynomial approximation of f(x): L ← p3. (l) If α3 > α, call the recursive procedure for the interval (x3,xb]. 3. Final step Add the last point pb to the polynomial approximation of f(x) : L ← pb and finish the combined sampling procedure. Geoinformatics FCE CTU 17(2), 2018 14 T. Bayer: Efficient plotting the functions with discontinuities Method M3. By summarizing the facts mentioned above the proposed algorithm uses a triple recursion. In each step, the refined f(x) is approximated by at least three new points. For the pseudocode, see Algs. 1, 2. Compared to [15], this solution has a dramatically im- proved performance and requires fewer subdivisions. Due to the significantly higher recursion depth d, for fast and efficient estimation of the f(x) curvature, a single recursive step is not sufficient. This issue is illustrated in Tab. 1; our method is labeled as M3, a single recursion step method as M1. 4. Singularity detection This section describes several simple strategies for handling and detecting the discontinu- ities; their overview can be found in [24], [2], [26]. Early methods are based on the Markov models [18], [28], [25]. Currently, there are several approaches: applications of the Fourier transformation [5], [12], [4], [23], [17], [16], [11], wavelets [35], [45], Chebyshev series [31], triangulations [19]. Geometric measures of the curvature are described in [38], [22], [34], [1], [30], [36], the statistical-based methods in [8], [26]. Unlike the solutions searching for all discontinuities of f at entire Ω, each sampled point xi is checked for a discontinuity; the removable, jump and infinite discontinuities are involved. These testing criteria are local, they analyze behavior of the function in a boundary B(xi,ε) of xi, covered by the equally spaced points xi−2 = xi −2h, xi−1 = xi −h, and, xi+1 = xi + h, xi+2 = xi + 2h, where h = ε/2. For practical computation ε = 0.001. An infinite discontinuity is detected if (|fi−k| > y) ∨ (fi−k ≡±Inf) ∨ (fi−k ≡ NaN) , k = −2, .., 2. where y is the given threshold, NaN and Inf are the symbols for the positive infinity and the result of the undefined operation. The remaining discontinuities may be found using the criteria measuring the smoothness by changes in the variation. Two criteria, WENO and LR, described in [30], are presented. The WENO criterion wj(x), j = 0, 1, 2, is written as wj(x) = αi∑2 j=0 αi , αj = 1 (ISj + �)2 , where IS0 = 13 12 (fi−2 − 2fi−1 + fi)2 + 1 4 (fi−2 − 4fi−1 + fi)2, IS1 = 13 12 (fi−1 − 2fi + fi+1)2 + 1 4 (fi−1 −fi+1)2, IS2 = 13 12 (fi+2 − 2fi+1 + fi)2 + 1 4 (fi+2 − 4fi+1 + fi)2. If w0(x) > 1/3∨w1(x) > 1/3∨w2(x) > 1/3, f is probably not smooth at xi. The LR criterion has the following form LR(x) = ∣∣f2r −f2l ∣∣ f2r + f2l , (4.1) where fr = 3fi −4fi+1 + fi+2, fl = 3fi −4fi−1 + fi−2. If LR > 0.8, f is probably not smooth at xi. For practical computations, the criteria provide similar results. However, the WENO criterion seems to be more sensitive, and a steep slope classifies as the jump discontinuity. This issue is widely discussed in [30]. Geoinformatics FCE CTU 17(2), 2018 15 T. Bayer: Efficient plotting the functions with discontinuities Table 1: The depth of recursion for different values of the refinement criteria α in the methods M1 and M3 (proposed). Method Refinement criterion α[ ◦] 20 15 10 5 2 1 0.5 0.2 0.1 M1 7 9 9 19 55 105 217 535 1071 M3 3 3 3 7 15 27 57 121 297 5. Combined sampling with the singularities The proposed method checks for a discontinuity at each sampled point xi ∈ Ω using the rules mentioned above, where the amount of discontinuities denoted as k is not a priori known. It does not represent a rigid mathematical solution describing the behavior of the analyzed function, but only a simple method applicable to the technical computing. Our approach follows a heuristic sufficient for most functions, especially for the meridian or parallel coordinate functions. As a result, the set of disjoint subsets Ωg, Ωg ⊆ Ω, containing “good” data that allows for adaptive sampling is constructed; this technique was used in [26]. The point xi is classified as “good” if no singularity at f(xi) occurs. An interval Ω containing only “good” points is classified as “good” and labeled as Ωg. Unfortunately, the procedure cannot be generalized for higher-dimensional problems. Suppose the j − th interval Ωj = [aj,bj] containing a singularity c, c ∈ Ωj, and ε, ε > 0, representing the numerical threshold. In general, several cases need to be distinguished: • Case 1: the coincidence with the lower bound If aj ≡ c, the discontinuity c coincides with the lower bound of the interval Ωj. • Case 2: the proximity to the lower bound If aj < c∧|c−aj| < ε, the discontinuity c is close the lower bound of the interval Ωj. • Case 3: the coincidence with the upper bound If bj ≡ c, the discontinuity c coincides with the upper bound of the interval Ωj. • Case 4: the proximity to the upper bound If bj > c∧|c− bj| < ε, the discontinuity c is close the upper bound of the interval Ωj. • Case 5: the interior singularity If aj < c < bj, where |c−aj| > ε∧|c− bj| > ε, the discontinuity c lies inside the interval Ωj. For practical computation in the floating-point arithmetic, Cases 1, 2 and Cases 3, 4 may be joined, and we search for the discontinuity close to the lower or upper bounds. The modified conditions are: • Cases 1+2: the proximity/coincidence to the lower bound If aj ≤ c∧|c−aj| < ε, the discontinuity c coincides or it is close to the lower bound of the interval Ωj. Geoinformatics FCE CTU 17(2), 2018 16 T. Bayer: Efficient plotting the functions with discontinuities • Cases 3+4: the proximity/coincidence to the upper bound If bj ≥ c∧|c− bj| < ε, the discontinuity c coincides or it is close to the upper bound of the interval Ωj. Another two essential cases need to be resolved: • Case 6: Too narrow interval If bj −aj < ε, an “empty” interval Ωj arises. • Case 7: The incorrect interval If aj > bj, the incorrect interval Ωj appears. In general, they occur due to the behavior of the sampled function as well as a result of the floating-point arithmetic. For our problem, the empty interval is “unpromising” and (probably) does not contain any important data1. Moreover, there is no chance of the possible improvement; the interval cannot be expanded. In most cases, the “incorrect” intervals are also empty, or they become empty during the next processing. In general, these intervals may be rejected from further processing; Cases 6, 7 can be solved simultaneously. Depending on the position of the singularity c, three types of operations are performed: • Delete empty/incorrect Ωj The empty or incorrect interval Ωj is deleted. • Shrinking Ωj The interval Ωj with the discontinuity c close to the bounds is shrunk from the left/right. For Cases 1+2, aj = aj + ε, for Cases 3+4, bj = bj −ε. • Splitting Ωj The interval Ωj with the internal discontinuity c is split so that the new disjoint intervals Ωj,1 = [aj,c−ε], and Ωj,2 = [c + ε,bj) are created. It is obvious that the Case 6 appears as the result of the incorrect split or shrink operations, while the Case 7 is the result of the incorrect shrink operation. 5.1. Combined sampling algorithm involving the singularities Let us summarize the facts mentioned above into the algorithm. Our implementation is based on the stack S, see Alg. 3, the amount of Ωj splits denoted s represents the recursion depth. The basic idea is to set Ωg ≡ Ω, to loop over all the adaptively sampled points pi, to check for a singularity c in the boundary B(xi,ε) of pi and to make a decision about the Ωj boundaries. If no singularity occurs, all points are classified as good, and the polygonal approximation Lj of the curve is constructed. All disjoint polygonal approximations Lj are stored in the list L. The discontinuities are localized successively, one by one. The stack-based approach consists of two phases (initial and recursive): 1However, for some specific types of non-continuous functions this data may be important (functions with isolated points) and cannot be skipped. Geoinformatics FCE CTU 17(2), 2018 17 T. Bayer: Efficient plotting the functions with discontinuities Algorithm 3 Combined sampling with the singularities, the stack implementation. 1: function csSingStack(f, L,aj,bj,s,d,d,ε,α) 2: S = ∅,s = 0 3: S ← Ωj = [aj,bj] 4: while S 6= ∅ do 5: Ωj = S.pop() 6: try 7: Lj = ∅ 8: csInit(f,Lj,aj,bj, 1, 1,d,d,ε,α) 9: L ← Lj 10: catch (SingularityException e) 11: c = e.x 12: if (s > s) then 13: L = ∅ 14: return 15: else 16: k = 0,Ωj,1 = Ωj,Ωj,2 = Ωj 17: processInt(Ωj,c,Ωj,1,Ωj,2,k,s,ε) 18: if (k > 0) then 19: S ← Ωj,1 20: else if (k > 1) then 21: S ← Ω2,1 1. The initial phase Initialize the empty stack S = ∅. Create Ωg = [a,b] and push S ← Ωg. 2. The recursive steps Repeat the following steps until S is empty: (a) Pop the actual good interval Ωgj ← S from S and get aj, bj. (b) Create the empty list Lj = ∅. (c) Create the temporary polygonal approximation of f on Ωj = [aj,bj] using com- bined sampling stored in Lj. (d) If no discontinuity appears, add Lj to L: L ← Lj and go to step (a). Otherwise, c represents the discontinuity that must be treated in Steps e-h). (e) If s > s, the maximum allowed recursion depth is exceeded without a reasonable solution. Clear the polygonal approximation L. (f) Otherwise, initialize the newly created intervals Ωj,1 = [aj,1,bj,1], Ωj,2 = [aj,2,bj,2], as Ωj,1 = Ωj, Ωj,2 = Ωj, and the amount of created intervals as k = 0. (g) Call the function processInt with parameters Ωj, c, Ωj,1, Ωj,2, k, S, ε refining the interval Ωj; see Alg. 4. The subintervals Ωj,1, Ωj,2 are passing by reference. Geoinformatics FCE CTU 17(2), 2018 18 T. Bayer: Efficient plotting the functions with discontinuities (h) If at least one new interval needs to be created, then k > 0. Push Ωj,1 to the stack: S ← Ωj,1. If k > 1, push the second interval Ωj,2 to the stack: S ← Ωj,2. Processing the interval. The procedure processInt(Ωj,c,Ωj,1,Ωj,2,k, S,ε) refines the interval Ωj. Depending on the c value, the Ωj bounds are shifted, or Ωj is split to Ωj,1, Ωj,2. It can be summarized as follows: 1. If aj > bj, Ωj represents an incorrect interval, return. 2. If |bj −aj| < ε, Ωj represents an “empty” interval, return. As mentioned above, for some functions with the multiple jump discontinuities the empty interval cannot be skipped. 3. If aj ≤ c∧|c−aj| < ε, the discontinuity c is close to the lower bound of the interval Ωj. Shift the lower bound aj,1 = aj + ε of Ωj,1 and increase the amount of created intervals k = k + 1. 4. If bj ≥ c∧|c− bj| < ε, the discontinuity c is close to the upper bound of the interval Ωj. Shift the upper bound bj,1 = bj −ε of Ωj,1 and increase the amount of created intervals k = k + 1. 5. If aj < c < bj, then |c−aj| > ε∧|c− bj| > ε, the discontinuity c is inside the interval Ωj which needs to be split to Ωj,1, Ωj,2. Shift the upper bound bj,1 = c − ε of Ωj,1 and the lower bound aj,2 = c + ε of Ωj,2. Increase the amount of the created intervals k = k + 2 and splits s = s + 1. 6. If at least one new interval needs to be created, then k > 0. Push Ωj,1 to the stack: S ← Ωj,1. If k > 1, push the second interval Ωj,2 to the stack: S ← Ωj,2. For the implementation, see Alg. 4. Analogously, the singularity value c is stored in the thrown exception, and processed in the try-catch block. In general, the stack-based imple- mentation is more stable than a common recursion and does not suffer from too high value of the recursion depth d; especially for the small values of ε. 5.2. Utilization in geoinformatics The proposed methods may be used for the polygonal approximation of curves when a more compact and efficient representation is required. The circles, circular arcs, ellipses or offsets of curves (e.g., buffers) cannot be internally stored in the *.shp files. It can also be used for the contour lines simplification (removing the adjacent segments, where αi < α). Another utilization for which the algorithm was originally developed, is a more efficient reconstruc- tion of the map projection graticule. The coordinate functions F(ϕ,λ), G(ϕ,λ) of the map projection depend on ϕ,λ and may contain several discontinuities. Therefore, the problem needs to be generalized, and its solution must be adapted for these facts. The algorithm for the map projection graticule construction will be presented in the next paper. 6. Experiments and results. The proposed methods for combined sampling have been tested for the set of 9 functions; the ability to detect and treat the discontinuities represents an important factor. Geoinformatics FCE CTU 17(2), 2018 19 T. Bayer: Efficient plotting the functions with discontinuities Algorithm 4 Processing the interval: shift bounds or split. 1: function processInt(Ωj,c,Ωj,1,Ωj,2,k,s,ε) 2: if (aj > bj) then 3: return 4: if |bj −aj| < ε then 5: return 6: if (aj ≤ c) ∧ (|c−aj| ≤ ε) then 7: aj,1 = aj + ε 8: k = k + 1 9: else if (bj ≥ c) ∧ (|c− bj| ≤ ε) then 10: bj,1 = bj −ε 11: k = k + 1 12: else if (aj < c < bj) ∧ (|c−aj| > ε) ∧ (|c− bj| > ε) then 13: bj,1 = c−ε,aj,2 = c + ε 14: k = k + 2 15: s = s + 1 6.1. List of functions In all cases is Ω = [−5π, 5π], ε = 0.001, α = 1◦, the thresholds s, d have not been set. The first function f1(x) = { 1 + x, x ≥ 0, 0, x < 0, has the jump discontinuity at x = 0. The polyline representation contains only 4 points divided into two intervals Ωg1 = [−5π,−1.3 · 10−4], Ω g 2 = [1.6 · 10−4, 5π], the amount of splits is s = 1, no recursion is required d = 0. The second function f2 (x) = e −500x2, is quite steep and does not contain any discontinuity. These facts led to the polygonal approximation formed by 266 points, Ωg = Ω, no splits, s = 0, the maximum recursion depth is d = 3. The third function f3 (x) = sin(x)/x, has a removable discontinuity at x = 0, where f(0) = 1. Ω is divided into two subintervals Ω g 1 = [−5π,−1.3 · 10−4], Ω g 2 = [1.6 · 10−4, 5π], the amount of splits is s = 1, no recursion required d = 0, the polygonal approximation contains 122 points. The fourth function f4(x) = x/ sin(x), has jump discontinuities at x = ±π ± kπ. Hence, Ω is split into 9 sub intervals: Ωg1 = [−15.692,−12.579], Ωg2 = [−12.554,−9.434], Ω g 3 = [−9.415,−6.290], Ω g 4 = [−6.277,−3.145], Ω g 5 = [−3.138, 3.138], Ω g 6 = [3.145, 6.277], Ω g 7 = [6.290, 9.415], Ω g 8 = [9.434, 12.554], Ω g 9 = [12.579, 15.692]. The amount of splits s = 8 as well as the recursion depth d = 6, provide the polygonal approximation containing 984 points. The next function f5(x) = x sin(5/x), Geoinformatics FCE CTU 17(2), 2018 20 T. Bayer: Efficient plotting the functions with discontinuities -10 0 10 x 0 5 10 15 20 f( x ) f1(x), 4 sampled points -10 0 10 x -10 -5 0 5 10 f( x ) f2(x), 266 sampled points -10 0 10 x -10 -5 0 5 10 f( x ) f3(x), 122 sampled points -10 0 10 x -10 -5 0 5 10 f( x ) f4(x), 984 sampled points -10 0 10 x -5 0 5 10 f( x ) f5(x), 6032 sampled points -10 0 10 x -10 -5 0 5 10 f( x ) f6(x), 755 sampled points -10 0 10 x -10 -5 0 5 10 f( x ) f7(x), 163 sampled points -10 0 10 x -10 -5 0 5 10 f( x ) f8(x), 1948 sampled points -10 0 10 x -10 -5 0 5 10 f( x ) f9(x), 1109 sampled points Figure 5.1: The polygonal approximation of functions f1(x), ...,f9(x) created by the combined sampling technique; the discontinuities involved. has a discontinuity at x = 0, and f(0) = 0. Ω should be divided in two subintervals, but the algorithm creates 12 subintervals in 42 splits: Ωg1 = [−15.708,−0.015], Ω g 2 = [−0.009,−0.008], Ω g 3 = [−0.007,−0.006], Ω g 4 = [−0.005,−0.004], Ω g 5 = [−0.004,−0.003], Ω g 6 = [−0.003, 0.003], Ω g 7 = [0.003, 0.004], Ω g 8 = [0.004, 0.005], Ω g 9 = [0.005, 0.006], Ω g 10 = [0.006, 0.007], Ω g 11 = [0.008, 0.009], Ωg12 = [0.015, 15.708]. As a result of the recursion depth d = 9, the total amount of points is n = 6032. While the graphical representation is aesthetically pleasing, the polygonal representation contains redundant data; this is due to the false detection of discontinuities provide by LR criterion. The next function f6(x) = 1/(tan 2x tan 0.5x), has the infinite singularities at x = ±1/2π ± kπ, and at x = 0 ± k2π (total 15). The adaptive sampling procedure with the recursion depth d = 7 creates 755 points in 16 new intervals Ωg1 = [−15.708,−14.138, Ω g 2 = [−14.137,−12.598], Ω g 3 = [−12.535,−10.996], Ω g 4 = [−10.995 − 7.855], Ωg5 = [−7.853,−6.315], Ω g 6 = [−6.251,−4.713], Ω g 7 = [−4.712,−1.571], Ω g 8 = [−1.570,−0.032], Ω g 9 = [0.032, 1.570], Ω g 10 = [1.571, 4.712], Ω g 11 = [4.713, 6.251], Ω g 12 = [6.315, 7.853], Ωg13 = [7.855, 10.995], Ω g 14 = [10.996, 12.535], Ω g 15 = [12.598, 14.137], Ω g 16 = [14.138, 15.708] with s = 15 splits. The next function f7(x) = e−2x/(x− 1), Geoinformatics FCE CTU 17(2), 2018 21 T. Bayer: Efficient plotting the functions with discontinuities Algorithm 5 Plotting the function f9(x) involving discontinuities in the Octave v. 4.4 scripting language. xmin = -5*pi; xmax = -xmin; eps = 0.001; lrmax = 0.8; x = xmin:0.001:xmax; f = @(x)exp(x)./tan(x); y = f(x); d = lr( f, x, eps) > lrmax; x(d==1) = nan; y(d==1) = nan; plot(x,y,’-r’); xlim([xmin,xmax]) ylim([xmin,xmax]) set(gca,’XLim’,[xmin xmax]) set(gca,’YLim’,[xmin xmax]]) xlabel(’x’); ylabel(’f(x)’); axis equal; has the infinite singularity at x = 1, Ω is divided into two subintervals Ωg1 = [−4.286, 1.000], Ω g 2 = [1.000, 15.708]. The amount of splits is s = 1 as well as the recursion depth d = 6, provide the polygonal approximation by 163 points. The next function f8(x) = x2 sin 2x/(2x− 1), has the infinity singularity at x = 0.5, Ω is divided into two subintervals Ωg1 = [−15.708, 0.500], Ω g 2 = [0.500, 15.708]. The amount of splits is s = 1 and the recursion depth d = 7 bring the polygonal approximation by 1948 points. The last function f9(x) = ex/ tan x, has infinite singularities at x = ±kπ (total 11). The adaptive sampling procedure with the recursion depth d = 6 creates 1109 points in 10 new intervals Ωg1 = [−15.708,−12.567], Ω g 2 = [−12.566,−9.425], Ωg3 = [−9.425,−6.283], Ω g 4 = [−6.2837,−3.142], Ω g 5 = [−3.141,−0.001], Ω g 6 = [0.001, 3.119], Ω g 7 = [3.165, 5.925], Ω g 8 = [7.224, 8.138], Ω g 9 = [10.979, 11.012], Ω g 10 = [14.137, 14.138], with s = 10 splits. The polygonal approximations of all analyzed functions can be found in Fig. 5.1. It is evident that their courses have been estimated correctly. However, for some functions, the proposed combined sampling algorithm brings less redundant representation (f5). In general, this method cannot compete with the well-known high-end solutions (Mathematica) involving the robust numerical techniques, but it may provide a reasonable representation of functions in technical computing. 6.2. Comparison with other systems For comparison, the functions will be plotted in two well-known systems. While the open- source software is represented by Octave (v. 4.4), the commercial software by Wolfram Mathematica (v. 11). Geoinformatics FCE CTU 17(2), 2018 22 T. Bayer: Efficient plotting the functions with discontinuities x f x x f x x f x x f x x f x x f x x f x x f x x f x f 1 (x) f 2 (x) f 3 (x) f 4 (x) f 5 (x) f 6 (x) f 7 (x) f 8 (x) f 9 (x) Figure 6.1: Functions f1(x), ...,f9(x) plotted in the Wolfram Mathematica software, v. 11. Octave. Unfortunately, the open-source software Octave (v. 4.4) does not support a correct plotting the functions involving discontinuities. The main idea of our solution is based on splitting Ω to the subintervals Ωg by putting NaN numbers between Ωg, where the function is undefined. The discontinuities are detected by the LR criterion given by Eq. 4.1. From a mathematic point of view, more characteristics of the function behavior need to be studied (first and second derivatives, asymptotes, local/global minima, ...), but this is outside the scope of the paper. The script can be found in Alg. 5. The previous results set the numerical characteristics: Ω = [−5π, 5π], ε = 0.001, LR = 0.8, sampling step h = 0.001. A curve is sampled uniformly by 10000π (approx. 31 000) points; this “relatively large” value has been set empirically. In general, the obtained results are analogous to the proposed method. Unfortunately, the algorithm is sensitive to the values of h, and LR. For the larger values of h, h = 0.01, only a subset of functions is sampled Geoinformatics FCE CTU 17(2), 2018 23 T. Bayer: Efficient plotting the functions with discontinuities Table 2: Uniform and combined sampling of the graticules of the equal area cylindrical, conic, azimuthal and Werner-Staab projections using combined sampling. The quantitative parame- ters are presented. Projection Sampling nmer npar αm αp α̃m α̃p Equal area cylindrical Uniform 1443 1406 0.00 0.00 0.00 0.00Combined 195 190 0.00 0.00 0.00 0.00 Equal area conic Uniform 1443 1406 0.00 2.20 0.00 2.20Combined 195 1612 0.00 5.01 0.00 1.89 Equal area azimuthal Uniform 1443 1406 0.00 5.00 0.00 5.00Combined 195 2476 0.00 4.74 0.00 2.81 Werner-Staab Uniform 1443 1406 9.27 5.00 3.86 2.93Combined 2334 1747 4.99 5.00 2.36 2.34 correctly. Mathematica. Wolfram Mathematica v. 11, the well-known system for technical comput- ing, developed for three decades, supports automatic plotting the functions with the disconti- nuities. The script has a straightforward form; see Alg. 6. For the functions f1(x), ...,f8(x) the analogous discontinuities have been detected. Unfortunately, for f9, the asymptote x = −3π has not been recognized; see Fig. 6.1. Moreover, the asymptote x = −2π is hard to distin- guish. If we understand Mathematica as the reference software, our proposed algorithm may be found as a simple and efficient tool for plotting a general function of the one variable. However, many situations may appear when it fails. Algorithm 6 Plotting the function f9(x) involving the discontinuities in the Mathematica v. 11 scripting language. xmin = -5*Pi; xmax = -xmin; F[x_] := Exp[x]/Tan[x]}; Plot[F[x], x, xmin, xmax, PlotRange -> xmin, xmax, xmin, xmax, AxesLabel -> x, f[x], AxesOrigin -> 0, 0, AspectRatio -> Automatic, PlotStyle -> Red]; 6.3. Construction of the map projection graticule In the last test, where the graticules of several map projections are reconstructed, the uniform and adaptive sampling techniques are compared regarding the data representation compact- ness. It is measured by the amount of the sampled meridian points nmer, and parallel points npar. Maximum angles αm,αp between the sampled meridian and parallel segments together with their mean values α̃m, α̃p are measured. The graticule is constructed over the entire planisphere, so Ω = Ωϕ×Ωλ, where ϕ ∈ [−π/2,π/2], and λ ∈ [−π,π], the offsets between the meridians and parallels are ∆ϕ = ∆λ = 10◦. In uniform sampling, the sampling steps of the meridians and parallels are δϕ = δλ = 2◦; combined sampling uses α = 2◦. For the circular arcs, uniform and combined sampling provide the analogous density of points representing the polygonal approximation. On the contrary, uniform sampling brings a higher density of the sampled points for the straight lines and vice versa for most of the curves. Four map Geoinformatics FCE CTU 17(2), 2018 24 T. Bayer: Efficient plotting the functions with discontinuities -18 0.0 -17 0.0 -16 0.0 -15 0.0 -14 0.0 -13 0.0 -12 0.0 -11 0.0 -10 0.0 -90 .0 -80 .0 -70 .0 -60 .0 -50 .0 -40 .0 -30 .0 -20 .0 -10 .0 0.0 10. 0 20. 0 30. 0 40. 0 50. 0 60. 0 70. 0 80. 0 90. 0 100 .0 110 .0 120 .0 130 .0 140 .0 150 .0 160 .0 170 .0 180 .0 -0.0 0.0 -90.0 -90.0 -80.0 -80.0 -70.0 -70.0 -60.0 -60.0 -50.0 -50.0 -40.0 -40.0 -30.0 -30.0 -20.0 -20.0 -10.0 -10.0 0.0 0.0 10.0 10.0 20.0 20.0 30.0 30.0 40.0 40.0 50.0 50.0 60.0 60.0 70.0 70.0 80.0 80.0 90.0 90.0 -180.0 -170.0 -160.0 -150.0 -140.0 -130.0 -120.0 -110.0 -100. 0 -90.0 -80.0 -70.0 -60. 0 -50. 0 -40. 0 -30. 0 -20. 0 -10. 0 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 110.0 120.0 130.0 140.0 150.0 160.0 170.0 180.0 -0.0 0.0 -90.0 -90.0 -80.0 -80.0 -70.0 -70.0 -60.0 -60.0 -50.0 -50 .0 -40.0 -40.0 -30.0 -30.0 -20.0 -20.0 -10.0 -10.0 0.0 0.0 10.0 10.0 20.0 20.0 30.0 30.0 40.0 40.0 50.0 50.0 60.0 60.0 70.0 70.0 80.0 80.0 90.0 90.0 -180.0 -170.0 -160.0 -150.0 -140.0 -130.0 -120.0 -110.0 -100.0 -90.0 -80.0 -70.0 -60.0 -50. 0 -40 .0 -30 .0 -20 .0 -10 .0 0.0 10 .0 20 .0 30. 0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 110.0 120.0 130.0 140.0 150.0 160.0 170.0 180.0 -0. 0 0.0 -90.0 -90 .0 -80.0 -80 .0 -70.0 -70 .0 -60.0 -60 .0 -50.0 -50 .0 -40.0 -40 .0 -30.0 -30 .0 -20.0 -20 .0 -10.0 -10 .0 0.0 0.0 10.0 10 .0 20.0 20 .0 30.0 30 .0 40.0 40 .0 50.0 50 .0 60.0 60 .0 70.0 70 .0 80.0 80 .0 90.0 90 .0 -180.0 -170.0 -160 .0 -150 .0 -140 .0 -13 0.0 -12 0.0 -11 0.0 -10 0.0 -90 .0 -80 .0 -70 .0 -60 .0 -50 .0 -40 .0 -30 .0 -20 .0 -10 .0 0.0 10 .0 20 .0 30 .0 40 .0 50 .0 60 .0 70 .0 80 .0 90 .0 10 0.0 110 .0 120 .0 130 .0 140 .0 150.0 160. 0 170.0 180.0 -0. 0 0.0 -90.0-90.0 -80.0 -80.0 -70.0 -70.0 -60.0 -60.0 -50.0 -50.0 -40.0 -40.0 -30.0 -30.0 -20.0 -20 .0 -10.0 -10. 0 0.0 0.0 10.0 10. 0 20.0 20 .0 30.0 30 .0 40.0 40 .0 50.0 50 .0 60.0 60 .0 70.0 70 .0 80.0 80 .0 90.0 90 .0 Figure 6.2: The reconstructed graticules of the equal area cylindrical, conic, azimuthal and Werner-Staab projections using the combined sampling technique. projections are involved in testing: equal-area cylindrical, conic (ϕ1 = 45◦), azimuthal, and Werner-Staab, their coordinate functions are continuous on Ω. The results are summarized in Tab. 2. For the straight segments, uniform sampling provides the redundant data; this issue refers to the cylindrical projection as well as to the meridians of the conic and azimuthal projections. While the constant curvature leads to the similar results (parallels of conic, azimuthal and Werner-Staab projections), in the higher-curvature regions, the combined sampling provides a smoother approximation (meridians of Werner-Staab projection). In general, combined sampling preserves the curvature better (see αm,αp, α̃m, α̃p values) and brings less redundant data which requires more sampled points. The reconstructed graticules can be found in Fig. 6.2. The modified version of the algorithm treating the discontinuities in the coordinate functions F,G will be presented in the next paper. 7. Conclusion This article presented a new algorithm combining the uniform and adaptive sampling tech- niques applicable to the functions involving the discontinuities, which are detected by the LR criterion. It can be used for the polygonal approximation of the curves when a more compact, efficient and less redundant representation is required. A typical example is represented by the circles, circular arcs, ellipses or offsets of curves, but it can be easily applicable to the map Geoinformatics FCE CTU 17(2), 2018 25 T. Bayer: Efficient plotting the functions with discontinuities projection graticule and similar problems in GIS and cartography. The illustrating examples indicate that the functions involving the discontinuities widely applied in technical practice can be plotted efficiently. However, there are many more complex functions, where the pro- posed solutions are not sufficient and need to be refined. A typical situation is represented by the isolated point, which is currently thrown. Another benefit is the relatively simple stack-based implementation. The source code written in Java can be found in the GitHub repository https://github.com/bayertom/sampling. The next paper brings an extension of this algorithm for combined sampling of the map projection graticules, when the coordinate functions F,G involve the discontinuities. References [1] Francesc Arandiga et al. “Interpolation and Approximation of Piecewise Smooth Func- tions”. In: SIAM Journal on Numerical Analysis 43.1 (2005), pp. 41–57. doi: 10.1137/ S0036142903426245. [2] Rick Archibald, Anne Gelb, and Jungho Yoon. “Determining the locations and disconti- nuities in the derivatives of functions”. In: Applied Numerical Mathematics 58.5 (2008), pp. 577–592. [3] André Ricardo Backes and Odemir Martinez Bruno. “Polygonal approximation of dig- ital planar curves through vertex betweenness”. In: Information Sciences 222 (2013), pp. 795–804. [4] Nana S. Banerjee, James F. Geer, and Langley Research Center. Exponential approxi- mations using fourier series partial sums [microform] / Nana S. Banerjee and James F. Geer. English. National Aeronautics and Space Administration, Langley Research Cen- ter ; National Technical Information Service, distributor Hampton, Va. : [Springfield, VA, 1997, p. 1 v. [5] Robert Bruce Bauer. “Band pass filters for determining shock locations”. In: (). [6] Frédéric Benhamou and William J. Older. “Applying interval arithmetic to real, integer, and boolean constraints”. In: The Journal of Logic Programming 32.1 (1997), pp. 1–24. issn: 0743-1066. doi: https://doi.org/10.1016/S0743-1066(96)00142-2. [7] P. Binev et al. “Adaptive approximation of curves”. In: Approximation Theory (2004), pp. 43–57. [8] Mira Bozzini and Milvia Rossini. “The detection and recovery of discontinuity curves from scattered data”. In: Journal of Computational and Applied Mathematics 240.Sup- plement C (2013). MATA 2012, pp. 148–162. issn: 0377-0427. doi: https://doi.org/ 10.1016/j.cam.2012.06.014. [9] A Carmona-Poyato et al. “Polygonal approximation of digital planar curves through break point suppression”. In: Pattern Recognition 43.1 (2010), pp. 14–25. [10] Filipe de Carvalho Nascimento et al. “Approximating implicit curves on plane and surface triangulations with affine arithmetic”. In: Computers & Graphics 40 (2014), pp. 36–48. Geoinformatics FCE CTU 17(2), 2018 26 https://github.com/bayertom/sampling https://doi.org/10.1137/S0036142903426245 https://doi.org/10.1137/S0036142903426245 https://doi.org/https://doi.org/10.1016/S0743-1066(96)00142-2 https://doi.org/https://doi.org/10.1016/j.cam.2012.06.014 https://doi.org/https://doi.org/10.1016/j.cam.2012.06.014 T. Bayer: Efficient plotting the functions with discontinuities [11] Dennis Cates and Anne Gelb. “Detecting derivative discontinuity locations in piecewise continuous functions from Fourier spectral data”. In: Numerical Algorithms 46.1 (2007), pp. 59–84. [12] Knut S Eckhoff. “Accurate reconstructions of functions of finite regularity from trun- cated Fourier series expansions”. In: Mathematics of Computation 64.210 (1995), pp. 671– 690. [13] M.H. van Emden. “Value Constraints in the CLP Scheme”. In: Constraints 2.2 (Oct. 1997), pp. 163–183. issn: 1572-9354. doi: 10.1023/A:1009705709733. [14] Richard Fateman. “Honest Plotting, Global Extrema, and Interval Arithmetic”. In: Papers from the International Symposium on Symbolic and Algebraic Computation. IS- SAC ’92. Berkeley, California, USA: ACM, 1992, pp. 216–223. isbn: 0-89791-489-9. doi: 10.1145/143242.143314. [15] Luiz Henrique de Figueiredo. “Adaptive sampling of parametric curves”. In: Graphics Gems V 5 (1995), pp. 173–178. [16] Anne Gelb and Eitan Tadmor. “Adaptive edge detectors for piecewise smooth data based on the minmod limiter”. In: Journal of Scientific Computing 28.2 (2006), pp. 279–306. [17] Anne Gelb and Eitan Tadmor. “Spectral reconstruction of piecewise smooth functions from their discrete data”. In: ESAIM: Mathematical Modelling and Numerical Analysis 36.2 (2002), pp. 155–175. [18] Stuart Geman and Donald Geman. “Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images”. In: IEEE Transactions on pattern analysis and machine intelligence 6 (1984), pp. 721–741. [19] Tim Gutzmer and Armin Iske. “Detection of discontinuities in scattered data approxi- mation”. In: Numerical Algorithms 16.2 (1997), pp. 155–170. [20] Timothy J. Hickey, Zhe Qju, and Maarten H. Van Emden. “Interval Constraint Plotting for Interactive Visual Exploration of Implicitly Defined Relations”. In: Reliable Com- puting 6.1 (Feb. 2000), pp. 81–92. issn: 1573-1340. doi: 10.1023/A:1009950630139. [21] Yifan Hu. “Efficient, high-quality force-directed graph drawing”. In: Mathematica Jour- nal 10.1 (2005), pp. 37–71. [22] Guang-Shan Jiang and Chi-Wang Shu. “Efficient implementation of weighted ENO schemes”. In: Journal of computational physics 126.1 (1996), pp. 202–228. [23] George Kvernadze. “Determination of the jumps of a bounded function by its Fourier series”. In: Journal of approximation theory 92.2 (1998), pp. 167–190. [24] David Lee. “Detection, Classification, and Measurement of Discontinuities”. In: SIAM Journal on Scientific and Statistical Computing 12.2 (1991), pp. 311–341. doi: 10 . 1137/0912018. [25] David Lee and Grzegorz W Wasilkowski. “Discontinuity detection and thresholding-a stochastic approach”. In: Journal of Complexity 9.1 (1993), pp. 76–96. [26] Licia Lenarduzzi and Robert Schaback. “Kernel-based adaptive approximation of func- tions with discontinuities”. In: Applied Mathematics and Computation 307.Supplement C (2017), pp. 113–123. issn: 0096-3003. doi: https://doi.org/10.1016/j.amc.2017. 02.043. Geoinformatics FCE CTU 17(2), 2018 27 https://doi.org/10.1023/A:1009705709733 https://doi.org/10.1145/143242.143314 https://doi.org/10.1023/A:1009950630139 https://doi.org/10.1137/0912018 https://doi.org/10.1137/0912018 https://doi.org/https://doi.org/10.1016/j.amc.2017.02.043 https://doi.org/https://doi.org/10.1016/j.amc.2017.02.043 T. Bayer: Efficient plotting the functions with discontinuities [27] Hélio Lopes, João Batista Oliveira, and Luiz Henrique de Figueiredo. “Robust adaptive polygonal approximation of implicit curves”. In: Computers & Graphics 26.6 (2002), pp. 841–852. [28] Jose Marroquin, Sanjoy Mitter, and Tomaso Poggio. “Probabilistic solution of ill-posed problems in computational vision”. In: Journal of the american statistical association 82.397 (1987), pp. 76–89. [29] Byron Nakos and V Miropoulos. “Local length ratio as a measure of critical points detection for line simplification”. In: Fifth Workshop on Progress in Automated Map Generalization, Paris, France. 2003. [30] M Oliveria et al. “Universal high order subroutine with new shock detector for shock boundary layer interaction”. In: In other words 10 (2009), pp. 1–2. [31] Ricardo Pachón, Rodrigo B Platte, and Lloyd N Trefethen. “Piecewise-smooth cheb- funs”. In: IMA journal of numerical analysis 30.4 (2009), pp. 898–916. [32] Afonso Paiva et al. “Approximating implicit curves on triangulations with affine arith- metic”. In: Graphics, Patterns and Images (SIBGRAPI), 2012 25th SIBGRAPI Con- ference on. IEEE. 2012, pp. 94–101. [33] Mohammad Tanvir Parvez and Sabri A Mahmoud. “Polygonal approximation of digital planar curves through adaptive optimizations”. In: Pattern Recognition Letters 31.13 (2010), pp. 1997–2005. [34] Leszek Plaskota and Grzegorz W Wasilkowski. “Adaption allows efficient integration of functions with unknown singularities”. In: Numerische Mathematik 102.1 (2005), pp. 123–144. [35] M. Rossini. “2D-discontinuity detection from scattered data”. In: Computing 61.3 (Sept. 1998), pp. 215–234. issn: 1436-5057. doi: 10.1007/BF02684351. [36] Yiqing Shen and Gecheng Zha. “Improvement of the WENO scheme smoothness esti- mator”. In: International Journal for Numerical Methods in Fluids 64.6 (2010), pp. 653– 675. issn: 1097-0363. doi: 10.1002/fld.2168. [37] Moshe Shpitalni, Yoram Koren, and CC Lo. “Realtime curve interpolators”. In: Computer- Aided Design 26.11 (1994), pp. 832–838. [38] Kaleem Siddiqi, Benjamin B Kimia, and Chi-Wang Shu. “Geometric shock-capturing ENO schemes for subpixel interpolation, computation, and curve evolution”. In: Com- puter Vision, 1995. Proceedings., International Symposium on. IEEE. 1995, pp. 437– 442. [39] C. Smith and N. Blachman. The Mathematica graphics guidebook. Addison-Wesley, 1995. isbn: 9780201826555. url: https : / / books . google . cz / books ? id = l % 5C _ MYAQAAIAAJ. [40] Jeffrey Allen Tupper. Graphing equations with generalized interval arithmetic. Univer- sity of Toronto, 1996. [41] Leland Wilkinson. “Algorithms for choosing the domain and range when plotting a function”. In: (1991), pp. 1–8. Geoinformatics FCE CTU 17(2), 2018 28 https://doi.org/10.1007/BF02684351 https://doi.org/10.1002/fld.2168 https://books.google.cz/books?id=l%5C_MYAQAAIAAJ https://books.google.cz/books?id=l%5C_MYAQAAIAAJ T. Bayer: Efficient plotting the functions with discontinuities [42] Daniel C.H. Yang and Tom Kong. “Parametric interpolator versus linear interpolator for precision CNC machining”. In: Computer-Aided Design 26.3 (1994). Special Issue:NC machining and cutter-path generation, pp. 225–234. issn: 0010-4485. doi: https : / / doi.org/10.1016/0010-4485(94)90045-0. [43] S-S Yeh and P-L Hsu. “The speed-controlled interpolator for machining parametric curves”. In: Computer-Aided Design 31.5 (1999), pp. 349–357. [44] Syh-Shiuh Yeh and Pau-Lo Hsu. “Adaptive-feedrate interpolation for parametric curves with a confined chord error”. In: Computer-aided design 34.3 (2002), pp. 229–237. [45] Tony Chan Zhou and HM Zhou. “Adaptive Eno-Wavelet Transforms For Discontinu- ous Functions”. In: CAM Report, No. 99-21, Dept. of Math., UCLA, Submit to SIAM Numer. Anal. Citeseer. 1999. Geoinformatics FCE CTU 17(2), 2018 29 https://doi.org/https://doi.org/10.1016/0010-4485(94)90045-0 https://doi.org/https://doi.org/10.1016/0010-4485(94)90045-0 Geoinformatics FCE CTU 17(2), 2018 30 T. Bayer: Efficient plotting the functions with discontinuities Introduction Related Work Combined sampling Polygonal approximation of the curve Combined sampling technique Singularity detection Combined sampling with the singularities Combined sampling algorithm involving the singularities Utilization in geoinformatics Experiments and results. List of functions Comparison with other systems Construction of the map projection graticule Conclusion