输入:离散的采样点坐标和高度值(x_0,y_0,value_0),(x_1,y_1,value_1)......(x_n, y_n, value_n)
输出:等高线图,如下所示
wiki上的Marching squares算法对此有很好的说明,我也是按照wiki上面的步骤来实现这个算法的,下面对该算法的步骤进行简要说明。
输入参数:
1、点的集合(x_0,y_0,value_0),(x_1,y_1,value_1)......(x_n, y_n, value_n) ;
2、高度值数组,即每条等高线表示的值,以上图为例是[0.0, 0.1, 0.2, .... 1.0]共10个区间,9个等高线值;
3、颜色数组,必须与高度值数组等长。
步骤:
1、构造矩阵
计算出输入点的x,y范围,即找到xmin, ymin, xmax, ymax,形成一个rectangle,然后拟定一个step,形成一个空的矩阵,step越小,等高线越平滑,计算速度也越慢。
2、矩阵插值
由输入点的值确定矩阵上每一个元素的值,这样一来,稀疏的输入点就变成了一个稠密的矩阵。一种简单的插值方式是:
matrix(x,y) = Sum(1/dist[i]^2 * value[i]) / Sum(1/dist[i]^2),其中dist[i]表示第i个点到(x,y)的欧式距离。
3、画等高线
到这里就可以完全按照wiki上的Marching squares来进行了。Marching squares算法把等高线绘制分成了两种情况,其一是只画线不上色,即isoLine算法,这种比较简单,不做讨论了;其二是画等高带,即isoBand算法,也就是本文采用的方法,如下说明:
对于每一个等高带,从高度数组中取出相临的两个阈值(如0.0, 0.1),从颜色数组中取一个颜色,用两个阈值把矩阵三值化成为一个只包含0,1,2的矩阵。0表示小于矩阵元素小于这两个阈值,1表示介于其中间,2表示大于这两个阈值。
这时,你的矩阵中就只有(0,1,2)三个值了,然后考虑矩阵中的每四个元素组成的单位正方形的状态,该状态由其四个角的元素值唯一确定,共3^4 = 81种,然后根据这81种情况分类讨论,分别上色了,此处请参考wiki的marching square页面。
如此,第一种颜色的等高带就画完了,把第3步的内容重复10次,画完所有的等高线。
补充一段第3步的actionscript代码,是我按照wiki的步骤实现的
![](https://i-blog.csdnimg.cn/blog_migrate/8f900a89c6347c561fdf2122f13be562.gif)
![](https://i-blog.csdnimg.cn/blog_migrate/961ddebeb323a10fe0623af514929fc1.gif)
public function isoband(data:qyMatrix, thArray:Array, dstLayer:GraphicsLayer, extent:Extent):Vector.<Polygon> { var x:int, y:int, k:int; var count:int=0; var squareWidth:Number = extent.width / (data.width - 1); var squareHeight:Number = extent.height / (data.height - 1); var dstPolygonVec:Vector.<Polygon> = new Vector.<Polygon>(thArray.length - 1); // 等值线的每一个阈值 for (k=0; k<thArray.length - 1; k++){ // 3值化 var stateMat:qyMatrix = new qyMatrix(data.width, data.height); for (y=0; y<data.height; y++){ for (x=0; x<data.width; x++){ if (data.getData(x, y) < thArray[k]){ stateMat.setData(x, y, 0); } else{ if (data.getData(x, y) < thArray[k+1]){ stateMat.setData(x, y, 1); } else{ stateMat.setData(x, y, 2); } } } } var x1:Number, y1:Number, x2:Number, y2:Number, x3:Number, y3:Number, x4:Number, y4:Number; var polygon:Polygon = new Polygon; for (y=0; y<data.height - 1; y++){ var ymin:Number = extent.ymin + y / (data.height - 1) * extent.height; var ymax:Number = extent.ymin + (y+1) / (data.height - 1) * extent.height; for (x=0; x<data.width - 1; x++){ var xmin:Number = extent.xmin + x / (data.width - 1) * extent.width; var xmax:Number = extent.xmin + (x+1) / (data.width - 1) * extent.width; // square 四角坐标 var p7:MapPoint = new MapPoint (xmin, ymin); var p9:MapPoint = new MapPoint (xmax, ymin); var p3:MapPoint = new MapPoint (xmax, ymax); var p1:MapPoint = new MapPoint (xmin, ymax); // square 四角数值 var d7:Number = data.getData(x, y); var d9:Number = data.getData(x+1, y); var d3:Number = data.getData(x+1, y+1); var d1:Number = data.getData(x, y+1); var mid:Number; // isoband的顶点坐标 var pt1:MapPoint = null; var pt2:MapPoint = null; var pt3:MapPoint = null; var pt4:MapPoint = null; var pt5:MapPoint = null; var pt6:MapPoint = null; var pt7:MapPoint = null; var pt8:MapPoint = null; var squareState:String = getSquareState(stateMat, x, y); switch (squareState) // total 81 cases { // no color case "2222": case "0000": break; // square case "1111": polygon.addRing([p7, p9, p3, p1]); break; // triangle 8 cases case "2221": x1 = p1.x + (thArray[k+1] - d1) / (d3 - d1) * squareWidth; y2 = p1.y - (thArray[k+1] - d1) / (d7 - d1) * squareHeight; pt1 = new MapPoint(x1, p1.y); pt2 = new MapPoint(p7.x, y2); polygon.addRing([pt1, p1, pt2]); break; case "2212": y1 = p3.y - (thArray[k+1] - d3) / (d9 - d3) * squareHeight; x2 = p3.x - (thArray[k+1] - d3) / (d1 - d3) * squareWidth; pt1 = new MapPoint(p9.x, y1); pt2 = new MapPoint(x2, p1.y); polygon.addRing([pt1, p3, pt2]); break; case "2122": x1 = p9.x - (thArray[k+1] - d9) / (d7 - d9) * squareWidth; y2 = p9.y + (thArray[k+1] - d9) / (d3 - d9) * squareHeight; pt1 = new MapPoint(x1, p7.y); pt2 = new MapPoint(p9.x, y2) polygon.addRing([pt1, p9, pt2]); break; case "1222": x1 = p7.x + (thArray[k+1] - d7) / (d9 - d7) * squareWidth; y2 = p7.y + (thArray[k+1] - d7) / (d1 - d7) * squareHeight; pt1 = new MapPoint(x1, p7.y); pt2 = new MapPoint(p7.x, y2); polygon.addRing([p7, pt1, pt2]); break; case "0001": x1 = p3.x - (thArray[k] - d3) / (d1 - d3) * squareWidth; y2 = p7.y + (thArray[k] - d7) / (d1 - d7) * squareHeight; pt1 = new MapPoint(x1, p1.y); pt2 = new MapPoint(p7.x, y2); polygon.addRing([pt1, p1, pt2]); break; case "0010": y1 = ymin + (thArray[k] - d9) / (d3 - d9) * squareHeight; x2 = xmin + (thArray[k] - d1) / (d3 - d1) * squareWidth; pt1 = new MapPoint(xmax, y1); pt2 = new MapPoint(x2, ymax); polygon.addRing([pt1, p3, pt2]); break; case "0100": x1 = p7.x + (thArray[k] - d7) / (d9 - d7) * squareWidth; y2 = p3.y - (thArray[k] - d3) / (d9 - d3) * squareHeight; pt1 = new MapPoint(x1, p9.y); pt2 = new MapPoint(p9.x, y2); polygon.addRing([pt1, p9, pt2]); break; case "1000": x1 = xmax - (thArray[k] - d9) / (d7 - d9) * squareWidth; y2 = ymax - (thArray[k] - d1) / (d7 - d1) * squareHeight; pt1 = new MapPoint(x1, ymin); pt2 =