Contour plotting for scattered data

Comput. & Graphics Vol. 11, No. 2, pp. 101-104, 1987

0097--8493/87 $3.00 + .00 © 1987 Pergamon Journals Ltd.

Printed in Great Britain.

Computer Graphics in India

CONTOUR PLOTTING FOR SCATTERED DATA D. G. SAWKAR,G. R. SHEVAREand S. P. KORUTHU • Aeronautical Engineering Department, Indian Institute of Technology, Bombay, India

INTRODUCTION

Contour plotting is one of the most important methods of data presentation. It is relatively simple for the data points arranged in a rectangular mesh. But for scattered data, contour plotting is not as simple a matter. In this case the scattered data are available in the form (x, y~, f,), i = 1 • • • N, with the only restriction that the points (x, Y3, i = 1 • • • N, are distinct. Various methods to plot contours for such data exist as of today [1 ]. The present work is an extension of the method developed by Shevare et aL [2], and is based on a method proposed by Neilson and Franke [3]. Neilson and Franke accomplish the task of contour plotting in three steps: (1) initial triangulation, (2) optimized triangulation, and (3) estimation of second derivatives in the chosen scheme of interpolation. Neilson and Franke [3] start contour plotting with an initial triangulation. They optimize it so as to produce a triangulation having minimum number of long thin triangles. This is required for certain error estimates which in turn involve interior angles of the triangles in the chosen scheme of interpolation. They use either max-min criterion or min-max criterion [4] for optimization of triangulation. Their study indicates that optimization affects only 5-10% of the sides of the triangles and max-min criterion gives approximately the same results as the rain-max criterion. Once optimized triangulation is ready they use various interpolation schemes with second order derivatives continuous for interpolation. THE

PROPOSED

SCHEME

In the present method the task of contour plotting is divided into four steps: (1) initial triangulation (Figs. l(a) & l(b)), (2) constrained optimized triangulation (Fig. l(c)), (3) subtriangulation (Fig. l(d)), and (4) interpolation. Each of these steps is described below in brief. (a) Initial triangulation There are many ways to triangulate the set of given data points. The method used here assumes the domain D as the convex hull of V~. Accordingly the convex hull is first obtained. The points lying on this convex hull are then excluded, and a new convex hull for the remaining points is obtained. This is continued until no more convex hulls can be drawn. The result of the procedure is shown in Fig. l(a). The points on adjacent convex hulls are then suitably joined [5] to obtain a mesh of triangles (Fig. l(b)).

(b) Constrained optimization Unlike in the case of optimization followed in [3], here we adopt a constrained optimization. In present scheme, triangles in two neighbouring convex hulls are considered at a time. The optimization procedure is carried out only for these triangles either by using maxmin criterion or rain-max criterion. After the optimization of the set of triangles within the two neighbouring convex hulls is completed, the set of triangles in the next two neighbouring convex hulls is considered. (c) Subtriangulation In the constrained optimized triangulation only a part of the total number of triangles are considered for optimization at a time; the procedure thus results in saving of computer time. However, it is obvious that the triangulation obtained is only near optimal. This disadvantage is completely offset by subtriangulation. Not only this, subtriangulation makes the task of interpolation much easier. In the suggested scheme, each triangle Tsjk is considered with one of its neighbours T,-jtas shown in Fig. 2(a). The vertices of these triangles are connected to produce a quadrilateral Q~m. Diagonal (i, j ) of Qokl is common to both the triangles. The value of the function at the point of intersection of two diagonals (i, j ) and (k, l) of Qukt is calculated by interpolating f i n Qok~.If the intersection of the diagonals lies between Vi and ~ then we have produced an additional pseudo data point. A triangle may have a maximum of 3 pseudo data points, thus giving rise to four triangles within the original triangle Tug. Alternatively, a triangle may not have any pseudo data point at all. Figs. 2(a)-2(d) show these cases of subtriangulation. A union of the set of original data points and the set of pseudo data points can be considered as a new set of given data points; the process of subtriangulation can thus be continued indefinitely. The calculation of the function at a pseudo data point needs some discussion. To start with let us consider a quadrilateral Q12.3,4. If the quadrilateral is square in shape there is an ambiguity as to which of the two possible diagonals one should ultimately choose. It should be realized that the ambiguity persists whether we use min-max criterion or max-min criterion. Subtriangulation helps us in resolving the ambiguity. In carrying out subtriangulation we assume that the function within the quadrilateral can be represented by a bilinear patch. This way we are using the entire information available with us for producing an additional 101

D. G. SAWKARet aL

102

(a) SUCCESSIVE

(bl CONEX

HULLS

IRIANGULATION(INITIAL)

o-

(c)

PSEUDO POINT

DAIA

(d)

CONS3"RAINE D OP 3"IMISE D "I'R IANGUL AI'ION

SUBTRIANGUL AT ION

Fig. 1. Triangulation. (pseudo) data point. Of course, as explained earlier question of generating the pseudo data point does not arise if the point of intersection of the diagonals falls out side the quadrilateral. Thus, it should be clear from this discussion that if subtriangulation is carried out it renders the process of optimized triangulation redundant. (d) Interpolation The interpolation of the function happens to be the simplest of the tasks. Let the value of the function for which the contour plot is desired be F0. The triangulation is searched till we get a triangle or a subtriangle Tijk with vertices V, Vj, Vk such that it contains the contour of value FO. Linear interpolation between the vertices of the triangle (or the subtriangle) gives two points which are the end points of the contour within the triangle (subtriangle). These two points are joined by a straight line. The process is continued till all the triangles have been checked. MERITS/DEMERITS OF THE METHOD

If z = f(x, y) happens to be to a smooth multiple valued function, the formation of successive convex hulls before the triangulations becomes essential. After the outermost convex hull is formed one can find two successive convex hulls within the outermost convex hull; one belonging to the region of the surface having negative z component of the outward unit normal whereas the other belonging to the region having pos-

itive z component of the outward unit normal. The successive convex hull pairs and the triangulation of the region between two successive convex hulls belonging to the same region of the surface can be carried out in parallel till the entire set of data points is exhausted. A scheme which straightaway proceeds with triangulation, i.e, not proceeding from the successive convex hulls to triangulation, may turn out to be efficient in cases but it will not be able to represent multiple valued functions for the obvious reasons. In our method we are advocating the use of only constrained optimisation. Otherwise for the sake of 5 10% of the triangles all the triangles need a search. And search is a costly procedure. In our scheme the number oftriangles needing search has been partitioned by the convex hulls, thereby making optimization faster. The number of triangles needing optimization after passing the test of constrained optimization could be well within 5%. The scheme of subtriangulation is a novel feature of the proposed method. Subtriangulation helps in producing smooth curves very effectively. The process produces a quad-tree structure similar to the classical Warnock's recursive algorithm. However, Warnock's recursive algorithm always produces four branches at every junction, whereas in our algorithm four or less than four branches may be available at every junction. Nonetheless, the intrinsic property of the algorithm reduces the computation time to a great extent; and the only triangles (or subtriangles) that contain the contour are required to be subdivided.

Contour plotting for scattered data

103

J

/ (a) PRODUCTION DATA

OF A

PSEUDO

(b)

POINT

( A Case of

A CASE

OF TWO

PSEUDO

DATA

POINTS single

pseudo

data

po int )

(c)

A CASE

OF T H R E E

PSEUDO

DATA POINTS.

(d)

A CASE

OF NO

PSEUDO

DATA

POINT.

Fig. 2. Generation of subtriangles. Diagonals producing pseudo data points (marked O) are shown by thin lines. Subtriangles are hatched. At every level of subtriangulation a diagonal is used to produce pseudo data point (and in turn the sides of subtriangles), if it lies within the quadrilateral. This renders the process of optimized triangulation unnecessary and the generated contours do not deviate from the accurate contours to any significant extent. Linear interpolation used in the present method offers three advantages: Firstly, it does not require a solution of simultaneous equations large in number; secondly, an error caused by a spurious data point is contained in only a small region; lastly, it must also be realized that getting the values of pairs (x, y) having the same value for z from a higher order interpolation requires a costly process of finding out the roots of either a quadratic or a cubic equation. RESULTS The developed method has been successfully used for producing the contour plots of a sphere. It turns out that the contours are the arcs of circles within a

close tolerance. The method has been used for producing the contours for the second data set used by Akima [6]. The produced contours reproduce the convexities and the concavities and match well with the results of [7]. CLOSURE The proposed scheme of contour plotting is preferable to many other existing methods (e.g., [3]) because of many reasons. 1. Our method uses only constrained optimized triangulation which is obviously faster than the process of optimized triangulation as the latter involves the scanning of all the triangles at the same time. 2. The proposed method uses only linear interpolation. This will result in saving of computer time required for the solution of a large number of linear algebraic equations with the second derivatives as the unknowns. It should be noted, however, that this does not result in the contours with broken lines as the

D. G. SAWKARet al.

i04

n u m b e r of stages in subtriangulation can be increased to suit the requirement of the desired smoothness. 3. The method suggested is fully "local"; thus, u n desirable undulations in the contour plots due to a spurious data point are contained in a small region. REFERENCES

1. Richard Franke, A critical comparison of some methods for interpolation of scattered data. Naval Postgraduate School, TR #NPS-53-79-003 (1979). 2. G. R. Shevare, D. G. Sawkar, and B. S. Babu, An interpolation method to scattered data. Presented in International Conference on Functional Analysis in Approximation Theory, Bombay, India (1985). 3. Gregory M. Neilson and Richard Franke, Surface con-

4. 5. 6. 7.

struction based upon triangulations. Surfaces in Computer Aided Design, Proceedings of a conference held at Mathematisches Forschung~nstitut Oberwnlfach F.G.R., April 25-30 (1982). I. Babuska and A. K. Aziz, On the angle condition in the finite element method, SIAM J. Numer. Anal, 13, pp. 214-226 (1976). B.A. Lewis and J. S. Robinson, Triangulation of planar regions with application. The Computer Journal 21, pp. 324-332 (1978). H. Akima, A method of bivariate interpolation and smooth surface fitting for irregularly distributed data points. ACM TOMS4, pp. 148-159 and 160-164 (1978). C. S. Petersen, Adaptive contouring of three-dimensional surfaces. Computer Aided Geometric Design l, pp. 6174 (1984).