# 3.1.2 使用绘图API绘制Contour的思路

A week or so back I wrote about a package I ported/modified to create the Delaunay triangulation in Flash with a few AS3 classes. As I noted there, such atriangulated irregular network (TIN) allows us to interpolateisolines — lines of constant value (aka isarithms, commonly called contours).

So, given a field of points (weather stations, say)…

…with one or more attributes attached (temperature, say)…

…a TIN can be constructed.

With the above TIN, values can be interpolated along each edge between the points of known values (control points). The interpolation is strictly linear (that is, the value 50 would be interpolated halfway along an edge whose control points were valued 48 and 52).

With a given contouring interval (I’m using 4 degrees F here), we can connect some of these interpolated points, creating our contour lines.

With the previous steps stripped away, this creates a passable isoline map.

The lines are rigid, though, and should be smoothed for presentation. I allow two methods for this. You can use the “simple” method, which just uses the built-in graphics method curveTo between the midpoint of each isoline segment (below with the isoline interval decreased to 3 degrees).

The above looks alright, but the curves are not continuous, closed loops can still have hard corners, and the isolines no longer pass through the interpolated points (we have therefore generalized an already-inaccurate interpolation). My compatriot Andy Woodruff, author of the glorious new Cartogrammar blog, offered to write a nice continuous curve method that ensured isolines would still pass through the interpolated values. You can read about the method inhis post. Here she blows:

Bringing it all together, then, and incorporating the only extra feature I wrote (tinting of isolines), here’s a nice finished isoline map of temperature across the U.S.

My new isolining package for Flash/ActionScript3 accomplishes all of the above, requiring only an array of point data with attribute values attached. The above example, was accomplished with the following lines of code (after drawing the U.S. states from a shapefile).

1. //first, generate the array of triangles (ITriangle objects) fromthe point data

2. var triangles:Array = Delaunay.triangulate(points);

3. Delaunay.drawDelaunay(triangles,points, triClip, false); //comment this out if youdon't want to draw the triangulation

4. //generate an array of isolines (isoline objects)

5. var isos:Array = IsoUtils.isoline(triangles, points, triClip,3, 0);

6. //create color and class arrays for tinting the isolines

7. var classesArray:Array = new Array(40, 44, 48, 52, 56, 60, 64, 68, 72, 76);

8. var colorsArray:Array = new Array(0x051CFD, 0x4602FD, 0x6D0EEB, 0x8400FF,0xC400FF, 0xEA00FF, 0xFF00E2, 0xFF0095, 0xFF0030, 0xFF0015, 0xFB3507);

9. //then, actually draw them, using a continuous curve

10. IsoUtils.drawIsolines(isos,triClip, "continuous",colorsArray, classesArray, .5, .95);

Keep in mind: triangulation is just one interpolation method, and is many ways the least technical (and accurate). More accurate interpolation techniques includeinverse-distance andkriging. ***If you’re having trouble, and your isoline interval is not an integer, check out the comment at line 171 ofisoUtils.as. Please fix that, BTW.

I meant to add other features, but since I started work this past week, I’m posting the package as-is, and invite others to modify. On my wishlist:

• hypsometric tinting, or color between the lines, would allow for more effective terrain or temperature mapping
• support for projections and other coordinate conversions in the drawIsolines method. I have packages for converting lat/long to a number of map projections, but currently the drawIsolines method doesn’t have support for passing a point coordinate conversion method.
• an animated demo. This thing’s lightning-fast, so why not?
• something that would be super wicked would be if someone would implementTanaka’s illuminated contours [pdf] method, that thickens/thins and darkens/lightens lines like so…

…creating beautiful relief maps like the one below

1）离散区域呈有限个三角单元；

2）绘制每个三角单元内的isoline；

3）遍历所有三角单元；

4）Done，Enjoy。

update: for a cool usage of Delaunay triangulation, see my isolining package for ActionScript 3

The Delaunay triangulation was invented in 1934 by Boris Delaunay. According to Paul Bourke,

The Delaunay triangulation is closely related geometrically to the Direchlet tesselation also known as the Voronoi or Theissen tesselations. These tesselations split the plane into a number of polygonal regions called tiles. Each tile has one sample point in its interior called a generating point. All other points inside the polygonal tile are closer to the generating point than to any other. The Delauney triangulation is created by connecting all generating points which share a common tile edge. Thus formed, the triangle edges are perpendicular bisectors of the tile edges.

（由于不会在CSDN插入flash，该Demo地在此：http://indiemaps.com/flash/clickTest_v1.swf，抄袭者注）

With the idea of creating client-side isolines and inspired in part by the spirit of Keith Peters’ ongoing MathWorld Problem of the Week, I set out to write an AS3 class to create the above triangulation. After an hour of hacking, I gave up, and decided to just port the algorithm from another language. I choseFlorian Jenett’s Java version, itself a port of Paul Bourke’s original C. The port was easy, and makes it a cinch to create a triangulation of, say, selected weather stations in the U.S.

```1. typedef struct {
2.  int p1,p2,p3;
3. } ITRIANGLE;
4. typedef struct {
5.  int p1,p2;
6. }IEDGE;
7. typedef struct {
8.  double x,y,z;
9. } XYZ;
10.
11. /*
12.    Triangulation subroutine
13.    Takes as input NV vertices in array pxyz
14.    Returned is a list of ntri triangular faces inthe array v
15.    These triangles are arranged in a consistentclockwise order.
16.    The triangle array 'v' should be malloced to 3 *nv
17.    The vertex array pxyz must be big enough to hold3 more points
18.    The vertex array must be sorted in increasing xvalues say
19.    qsort(p,nv,sizeof(XYZ),XYZCompare);
20.       :
21.    int XYZCompare(void *v1,void *v2)
22.    {
23.       XYZ *p1,*p2;
24.       p1 = v1;
25.       p2 = v2;
26.       if (p1->x < p2->x)
27.          return(-1);
28.       else if (p1->x >p2->x)
29.          return(1);
30.       else
31.          return(0);
32.    }
33. */
34. int Triangulate(int nv,XYZ *pxyz,ITRIANGLE *v,int *ntri)
35. {
36.  int *complete = NULL;
37.    IEDGE*edges = NULL;
38.  int nedge = 0;
39.  int trimax,emax = 200;
40.  int status = 0;
41.
42.  int inside;
43.  int i,j,k;
44.  double xp,yp,x1,y1,x2,y2,x3,y3,xc,yc,r;
45.  double xmin,xmax,ymin,ymax,xmid,ymid;
46.  double dx,dy,dmax;
47.
48.  /* Allocate memory for the completeness list, flag for eachtriangle */
49.    trimax= 4 * nv;
50.  if ((complete = malloc(trimax*sizeof(int))) == NULL) {
51.       status= 1;
52.  goto skip;
53.    }
54.
55.  /* Allocate memory for the edge list */
56.  if ((edges = malloc(emax*(long)sizeof(EDGE))) == NULL) {
57.       status= 2;
58.  goto skip;
59.    }
60.
61.  /*
62.       Find the maximum and minimumvertex bounds.
63.       This is to allow calculationof the bounding triangle
64.    */
65.    xmin= pxyz[0].x;
66.    ymin= pxyz[0].y;
67.    xmax= xmin;
68.    ymax= ymin;
69.  for (i=1;i<nv;i++) {
70.  if (pxyz[i].x < xmin) xmin = pxyz[i].x;
71.  if (pxyz[i].x > xmax) xmax = pxyz[i].x;
72.  if (pxyz[i].y < ymin) ymin = pxyz[i].y;
73.  if (pxyz[i].y > ymax) ymax = pxyz[i].y;
74.    }
75.    dx= xmax - xmin;
76.    dy= ymax - ymin;
77.    dmax= (dx > dy) ? dx : dy;
78.    xmid= (xmax + xmin) / 2.0;
79.    ymid= (ymax + ymin) / 2.0;
80.
81.  /*
82.       Set up the supertriangle
83.       This is a triangle whichencompasses all the sample points.
84.       The supertriangle coordinatesare added to the end of the
85.       vertex list. The supertriangleis the first triangle in
86.       the triangle list.
87.    */
88.    pxyz[nv+0].x= xmid - 20 * dmax;
89.    pxyz[nv+0].y= ymid - dmax;
90.    pxyz[nv+0].z= 0.0;
91.    pxyz[nv+1].x= xmid;
92.    pxyz[nv+1].y= ymid + 20 * dmax;
93.    pxyz[nv+1].z= 0.0;
94.    pxyz[nv+2].x= xmid + 20 * dmax;
95.    pxyz[nv+2].y= ymid - dmax;
96.    pxyz[nv+2].z= 0.0;
97.    v[0].p1= nv;
98.    v[0].p2= nv+1;
99.    v[0].p3= nv+2;
100.   complete[0] = FALSE;
101.   *ntri = 1;
102.
103. /*
104.      Includeeach point one at a time into the existing mesh
105.   */
106. for(i=0;i<nv;i++) {
107.
108.      xp =pxyz[i].x;
109.      yp =pxyz[i].y;
110.      nedge =0;
111.
112. /*
113.         Setup the edge buffer.
114.         Ifthe point (xp,yp) lies inside the circumcircle then the
115.         threeedges of that triangle are added to the edge buffer
116.         andthat triangle is removed.
117.      */
118. for (j=0;j<(*ntri);j++) {
119. if (complete[j])
120. continue;
121.         x1= pxyz[v[j].p1].x;
122.         y1= pxyz[v[j].p1].y;
123.         x2= pxyz[v[j].p2].x;
124.         y2= pxyz[v[j].p2].y;
125.         x3= pxyz[v[j].p3].x;
126.         y3= pxyz[v[j].p3].y;
127.         inside= CircumCircle(xp,yp,x1,y1,x2,y2,x3,y3,&xc,&yc,&r);
128. if (xc < xp && ((xp-xc)*(xp-xc)) > r)
129.            complete[j]= TRUE;
130. if (inside) {
131. /* Check that we haven't exceeded the edge list size */
132. if (nedge+3 >= emax) {
133.               emax+= 100;
134. if ((edges = realloc(edges,emax*(long)sizeof(EDGE)))== NULL) {
135.                  status= 3;
136. goto skip;
137.               }
138.            }
139.            edges[nedge+0].p1= v[j].p1;
140.            edges[nedge+0].p2= v[j].p2;
141.            edges[nedge+1].p1= v[j].p2;
142.            edges[nedge+1].p2= v[j].p3;
143.            edges[nedge+2].p1= v[j].p3;
144.            edges[nedge+2].p2= v[j].p1;
145.            nedge+= 3;
146.            v[j]= v[(*ntri)-1];
147.            complete[j]= complete[(*ntri)-1];
148.            (*ntri)--;
149.            j--;
150.         }
151.      }
152.
153. /*
154.         Tagmultiple edges
155.         Note:if all triangles are specified anticlockwise then all
156.               interioredges are opposite pointing in direction.
157.      */
158. for (j=0;j<nedge-1;j++) {
159. for (k=j+1;k<nedge;k++) {
160. if ((edges[j].p1 == edges[k].p2) &&(edges[j].p2 == edges[k].p1)) {
161.               edges[j].p1= -1;
162.               edges[j].p2= -1;
163.               edges[k].p1= -1;
164.               edges[k].p2= -1;
165.            }
166. /* Shouldn't need the following, see note above */
167. if ((edges[j].p1 == edges[k].p1) &&(edges[j].p2 == edges[k].p2)) {
168.               edges[j].p1= -1;
169.               edges[j].p2= -1;
170.               edges[k].p1= -1;
171.               edges[k].p2= -1;
172.            }
173.         }
174.      }
175.
176. /*
177.         Formnew triangles for the current point
178.         Skippingover any tagged edges.
179.         Alledges are arranged in clockwise order.
180.      */
181. for (j=0;j<nedge;j++) {
182. if (edges[j].p1 < 0 || edges[j].p2 < 0)
183. continue;
184. if ((*ntri) >= trimax) {
185.            status= 4;
186. goto skip;
187.         }
188.         v[*ntri].p1= edges[j].p1;
189.         v[*ntri].p2= edges[j].p2;
190.         v[*ntri].p3= i;
191.         complete[*ntri]= FALSE;
192.         (*ntri)++;
193.      }
194.   }
195.
196. /*
197.      Removetriangles with supertriangle vertices
198.      Theseare triangles which have a vertex number greater than nv
199.   */
200. for(i=0;i<(*ntri);i++) {
201. if (v[i].p1 >= nv || v[i].p2 >= nv || v[i].p3>= nv) {
202.         v[i]= v[(*ntri)-1];
203.         (*ntri)--;
204.         i--;
205.      }
206.   }
207.
208.skip:
209.   free(edges);
210.   free(complete);
211. return(status);
212.}
213.
214./*
215.   Return TRUE ifa point (xp,yp) is inside the circumcircle made up
216.   of the points(x1,y1), (x2,y2), (x3,y3)
217.   Thecircumcircle centre is returned in (xc,yc) and the radius r
218.   NOTE: A pointon the edge is inside the circumcircle
219.*/
220.int CircumCircle(doublexp,double yp,
221. doublex1,double y1,doublex2,double y2,doublex3,double y3,
222. double*xc,double *yc,double*rsqr)
223.{
224. doublem1,m2,mx1,mx2,my1,my2;
225. doubledx,dy,drsqr;
226. doublefabsy1y2 = fabs(y1-y2);
227. doublefabsy2y3 = fabs(y2-y3);
228.
229. /*Check for coincident points */
230. if(fabsy1y2 < EPSILON && fabsy2y3 < EPSILON)
231. return(FALSE);
232.
233. if(fabsy1y2 < EPSILON) {
234.      m2 = -(x3-x2) / (y3-y2);
235.      mx2 =(x2 + x3) / 2.0;
236.      my2 =(y2 + y3) / 2.0;
237.      *xc =(x2 + x1) / 2.0;
238.      *yc = m2* (*xc - mx2) + my2;
239.   } elseif (fabsy2y3 < EPSILON) {
240.      m1 = -(x2-x1) / (y2-y1);
241.      mx1 =(x1 + x2) / 2.0;
242.      my1 =(y1 + y2) / 2.0;
243.      *xc =(x3 + x2) / 2.0;
244.      *yc = m1* (*xc - mx1) + my1;
245.   } else{
246.      m1 = -(x2-x1) / (y2-y1);
247.      m2 = -(x3-x2) / (y3-y2);
248.      mx1 =(x1 + x2) / 2.0;
249.      mx2 =(x2 + x3) / 2.0;
250.      my1 =(y1 + y2) / 2.0;
251.      my2 =(y2 + y3) / 2.0;
252.      *xc =(m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
253. if (fabsy1y2 > fabsy2y3) {
254.         *yc= m1 * (*xc - mx1) + my1;
255.      } else {
256.         *yc= m2 * (*xc - mx2) + my2;
257.      }
258.   }
259.
260.   dx = x2 - *xc;
261.   dy = y2 - *yc;
262.   *rsqr = dx*dx + dy*dy;
263.
264.   dx = xp - *xc;
265.   dy = yp - *yc;
266.   drsqr = dx*dx + dy*dy;
267.
268. //Original
269. //return((drsqr<= *rsqr) ? TRUE : FALSE);
270. //Proposed by Chuck Morris
271. return((drsqr- *rsqr) <= EPSILON ? TRUE : FALSE);
272.}```

0 条评论

• ### 4.1 数值积分、高等函数绘制

is defined informally as the signed area of the region in the xy-plane that is ...

• ### 3.1.3 绘制三维Contour图的思路

2007年秋，开始接触数值计算，看到Contour图形，我觉得很神奇，很好看。强烈的好奇心驱使下，零零碎碎看了相关文献，都看不懂。大约2009年深秋，我读到的最...

• ### 天太热画一个泵与风机给自己降温

按了水泵电源按钮后，上水阀门旋转开启，水流入水池，流入到一定高度时，管道里开始有水流动，测压管水位缓缓上升，集水池水位缓缓增加，最后回流到潜水泵箱体。关...

• ### 从MapX到MapXtreme2004[3]－搜索图元Feature

一、根据名称搜索图元 　　1、Mapxtreme的架构和Mapx有所变化，Mapx中，Layer包含Features，而Mapxtreme中则不是 　　2、...

• ### Codeforces Beta Round #51 C. Pie or die（博弈 思维）

Volodya and Vlad play the following game. There are k pies at the cells of n  × ...

• ### PAT (Advanced Level) Practice 1099 Build A Binary Search Tree （30 分）

A Binary Search Tree (BST) is recursively defined as a binary tree which has the...

• ### 如何使用CRM Document Template Profile

版权声明：本文为博主原创文章，遵循 CC 4.0 by-sa 版权协议，转载请附上原文出处链接和本声明。

• ### Spring_总结_04_高级配置(二)之条件注解@Conditional

在上一节，我们了解到 Profile 为不同环境下使用不同的配置提供了支持，那么Profile到底是如何实现的呢？其实Profile正是通过条件注解来实现的。

• ### 计算机如何理解我们的语言？NLP is fun！

【导读】我们从日常每天都会用到的推荐系统到现在研究火热的开放性聊天、对话机器人，越来越多的产品与应用的背后都需要自然语言处理（NLP）和知识图谱的技术。也有越来...