Efficient plotting the functions with discontinuities based on combined sampling

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. During the combined sampling, the polygonal approximation with a given threshold α between the adjacent segments is constructed. The automated detection and treatment of the discontinuities 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.


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 Ω g k without the internal singularities, 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 proposed 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 Ω g k , 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

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.

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 interpolation 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 polygonal 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.

Combined sampling
In this section, the proposed combined sampling technique providing the polygonal approximation 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.

Polygonal approximation of the curve
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 x i with the step δ, where δ = x i+1 −x i , provides a good approximation only if δ → 0. For the straight parts of the curve, many Geoinformatics FCE CTU 17(2), 2018 T. Bayer: Efficient plotting the functions with discontinuities US: α =10°, 63 sampled points 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.

Combined sampling technique
By avoiding colinear segments as well as a better adaptation to the different curvature, adaptive 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 illustrated in Fig. 3.1.Uniform sampling requires more data to maintain the same curvature represented by the angle α i between the adjacent segments (p i−1 , p i ) and (p i , p i+1 ).The difference 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 T. Bayer: Efficient plotting the functions with discontinuities 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 (p i−1 , p i , p i+1 ) to be three consecutive sampled points of the curve.The primary criterion is represented by the angular difference of both segments Different criteria can provide the analogous results: recall the distance of p i from p i−1 , p i+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.
if discontinuity in a then

The recursive step
Enter the recursive procedure and do the following substeps: into the approximate quarters where r 1 , r 2 , r 3 are the random numbers inside the interval [0.45, 0.55].This step prevents a situation, when ), but their intermediate 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 x 1 , x 2 , or x 3 occurs, throw the new exception with the argument indicating the singularity.

T. Bayer: Efficient plotting the functions with discontinuities
Algorithm 2 Combined sampling, the recursive phase.

Final step
Add the last point p b to the polynomial approximation of f (x) : L ← p b and finish the combined sampling procedure.
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 improved 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.
Unlike the solutions searching for all discontinuities of f at entire Ω, each sampled point x i 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(x i , ε) of x i , covered by the equally spaced points , where h = ε/2.For practical computation ε = 0.001.An infinite discontinuity is detected if 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 w j (x), j = 0, 1, 2, is written as where The LR criterion has the following form where If LR > 0.8, f is probably not smooth at x i .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].

Combined sampling with the singularities
The proposed method checks for a discontinuity at each sampled point x i ∈ Ω 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 x i is classified as "good" if no singularity at f (x i ) 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 = [a j , b j ] containing a singularity c, c ∈ Ω j , and ε, ε > 0, representing the numerical threshold.In general, several cases need to be distinguished: • Another two essential cases need to be resolved: • Case 6: Too narrow interval If b j − a j < ε, an "empty" interval Ω j arises.

• Case 7: The incorrect interval
If a j > b j , 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, a j = a j + ε, for Cases 3+4, b j = b j − ε.
• Splitting Ω j The interval Ω j with the internal discontinuity c is split so that the new disjoint intervals and 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.

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 p i , to check for a singularity c in the boundary B(x i , ε) of p i and to make a decision about the Ω j boundaries.If no singularity occurs, all points are classified as good, and the polygonal approximation L j of the curve is constructed.All disjoint polygonal approximations L j 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): T. Bayer: Efficient plotting the functions with discontinuities Algorithm 3 Combined sampling with the singularities, the stack implementation.

The initial phase
Initialize the empty stack S = ∅.Create Ω g = [a, b] and push S ← Ω g .

The recursive steps
Repeat the following steps until S is empty: (a) Pop the actual good interval Ω g j ← S from S and get a j , b j .
(b) Create the empty list L j = ∅.
(c) Create the temporary polygonal approximation of f on Ω j = [a j , b j ] using combined sampling stored in L j .
(d) If no discontinuity appears, add L j to L: L ← L j 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.
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 .
It can be summarized as follows: 1.If a j > b j , Ω j represents an incorrect interval, return.
2. If |b j − a j | < ε, Ω 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 a j ≤ c ∧ |c − a j | < ε, the discontinuity c is close to the lower bound of the interval Ω j .Shift the lower bound a j,1 = a j + ε of Ω j,1 and increase the amount of created intervals and increase the amount of created intervals and the lower bound a j,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 implementation 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 ε.

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 reconstruction 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.

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.
T. Bayer: Efficient plotting the functions with discontinuities 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: Ω g T. Bayer: Efficient plotting the functions with discontinuities

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 T. Bayer: Efficient plotting the functions with discontinuities Mathematica.Wolfram Mathematica v. 11, the well-known system for technical computing, developed for three decades, supports automatic plotting the functions with the discontinuities.The script has a straightforward form; see Alg. 6.For the functions f 1 (x), ..., f 8 (x) the analogous discontinuities have been detected.Unfortunately, for f 9 , the asymptote x = −3π has not been recognized; see Fig. 6.1.Moreover, the asymptote x = −2π is hard to distinguish.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.AxesOrigin -> 0, 0, AspectRatio -> Automatic, PlotStyle -> Red];

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 compactness.It is measured by the amount of the sampled meridian points n mer , and parallel points n par .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 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.

Conclusion
This article presented a new algorithm combining the uniform and adaptive sampling techniques 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 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 proposed 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.
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

Suppose d to beAlgorithm 1 1 :
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 T. Bayer: Efficient plotting the functions with discontinuities Combined sampling, the initial phase.function csInit(f, L, a, b, d, d, d, , α)

L 9 : 1 .
← P oint(a, y a ) as(f, L, a, b, y a , y b , d, d, d, , α) 10: L ← P oint(b, y b ) each pair of four new consecutive segments, the refinement criterion α k , k = 1, ..., 3, is evaluated 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 = {p i } n i=1 be the polynomial approximation of f (x) and Ω = [a, b] a subdomain.The algorithm may be summarized as follows: The initial phase Let L = ∅ be the empty set.Compute y a = f (a) and y b = f (b).If a singularity in a or b is detected, throw the exception.Add the initial vertex to L: L ← p a , where p a = (x a , y a ).Set the recursion depth d = 1.
(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

then 21 :
cs(f, L, x 3 , x b , y 3 , y b , d + 1, d, d, α); (e) Check the refinement criteria α 1 = α(p a , p 1 , p 2 ), α 2 = α(p 1 , p 2 , p 3 ), α 3 = α(p 2 , p 3 , p b ), 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, x 1 ).(g) Add new point p 1 to the polynomial approximation of f (x): L ← p 1 .(h) If α 1 > α ∨ α 2 > α, call the recursive procedure with the increased depth d = d + 1 for the interval (x 1 , x 2 ).(i) Add new point p 2 to the polynomial approximation of f (x): L ← p 2 .(j) If α 2 > α ∨ α 3 > α, call the recursive procedure with the increased depth d = d + 1 for the interval (x 2 , x 3 ).(k) Add new point p 3 to the polynomial approximation of f (x): L ← p 3 .(l) If α 3 > α, call the recursive procedure for the interval (x 3 , x b ].

Case 1 :
the coincidence with the lower bound If a j ≡ c, the discontinuity c coincides with the lower bound of the interval Ω j .• Case 2: the proximity to the lower bound If a j < c ∧ |c − a j | < ε, the discontinuity c is close the lower bound of the interval Ω j .• Case 3: the coincidence with the upper bound If b j ≡ c, the discontinuity c coincides with the upper bound of the interval Ω j .• Case 4: the proximity to the upper bound If b j > c ∧ |c − b j | < ε, the discontinuity c is close the upper bound of the interval Ω j .• Case 5: the interior singularity If a j < c < b j , where |c − a j | > ε∧|c − b j | > ε, 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 a j ≤ c ∧ |c − a j | < ε, the discontinuity c coincides or it is close to the lower bound of the interval Ω j .T. Bayer: Efficient plotting the functions with discontinuities • Cases 3+4: the proximity/coincidence to the upper bound If b j ≥ c ∧ |c − b j | < ε, the discontinuity c coincides or it is close to the upper bound of the interval Ω j .

Figure 6 . 2 :
Figure 6.2:The reconstructed graticules of the equal area cylindrical, conic, azimuthal and Werner-Staab projections using the combined sampling technique.

Table 1 :
The depth of recursion for different values of the refinement criteria α in the methods M1 and M3 (proposed).

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 parameters are presented.