3.1.2 使用绘图API绘制Contour的思路
3.1.2 使用绘图API绘制Contour的思路
原文：isolining package for ActionScript 3
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
If you add anything to the package, feel free to post a link to your revised version in the comments.
虽然没有细看算法，但是看来该大师做法也是：
1）离散区域呈有限个三角单元；
2）绘制每个三角单元内的isoline；
3）遍历所有三角单元；
4）Done，Enjoy。
但是他的三角单元如何离散得到，且看大师另一篇博文：delaunay triangulation in ActionScript 3，讲述了delaunay算法的前世今生，照抄如下：
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.
以下是C语言版本的该算法：
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.}本文分享自 传输过程数值模拟学习笔记 微信公众号，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ，欢迎热爱写作的你一起参与！
腾讯云开发者
