等值线的追踪算法

   在等值线的追踪中有一类相对来说比较简单的，也就是涉及不是很深的问题，那就是本文要说的基于规则的等值线的追踪网格。
　　先说一下等值线的形成，等值线的绘制是有两种方法的，一种就是完全测量，也就是在实际的施工中找到所有的等值点，标记上x,y,value，还有一种就是按一定的方式用几个预测点来采用一些插值方法形成规则的网格，每一个网格点上都有着坐标和value,　也就是高程值，然后在用这些网格上的值来估计出等值点，显然用第一种方法是不好的，常常也是不切实际的；这样也就有了一系列的等值线的追踪方法的形成；在这里，也有一个不得不说的就是基于几个观测点的网格化问题，这里面也有很多的值得研究的地方，
　　我在做项目时用的是Kriging插值算法，当然还有一些其它的很多，有加权反距离，有最小曲率等等很多，这些算法都是比较不好做的，
　　大家也可以到网上搜的了解一下，我自己也没有弄好它们，也是惭愧，在项目中用到的Kriging算法也是别人老师做的，但也没有做的非常的好，这个Kriging算法也有做的比较好的，我们做的是地质软件，象我们参照过的Srufer，就是这方面比较做的比较好 的，地质上也是用的比较多的，它的这方面的算法就做的比较好，可是它也还有一点问题，也就是我前面说过的，它也没有把断层做好，当然我也没有说别人的资格，我其实还远远没有做到他们的软件那样；在Kriging算法方面，也有一些开源的，比较好的就是斯坦福大学做的开源的Fortran的代码，那可真是利害，也做了开源软件的最高境界， 我开源了，你看不懂，等你看懂了，也过时了；他做的就是我开源了，你看不懂，那个可是真的太难看了，顺便说一句，如果大家有兴趣那个也真是值得大家研究研究，研究好了，那起码硕士论文那绝对是不成问题的，因为别人的这方面的博士论文也就那样，呵呵。　　以下就开始等值线追踪算法的讲解；等值线可以分成为非封闭等值线和封闭等值线，等值线的追踪先从网格边界或者网格内部上一等值点出发， 求得下一等值点，然后以此点出发搜索下一等值点（主要是求出坐标），一直这样下去，若遇到网格边界或者又回到起点（封闭等值线），则说明已经找到了一条等值上点所在的位置， 也就是找到了等值线。
　　要注意网格，什么是网格？网格又意味着什么？等值线在网格中又有什么不一样的地方？大家可要知道，我现在要说的是一种基于规则网格的等值线的追踪，就是那种矩形的网格，只有网格上才有坐标和高程值，而在算法中又要找出坐标，难道是找那些网格点，把网格点上的等值点长出来，这当然不是，我们要找的点应该是在那些网格的每一个小格子的边上，在遇到了等值点正好在网格点上时我们还不得不做一些小的处理，否则是无法追踪下去的，待会我再讲原因； 大家又可能会问，那没有值我们追踪个啥，
　　这也就是这种方法要做的，记住在每一个小格上的四个点的坐标和高程值我们是知道的，那就可以根据这些值来估计了，有人会说，是估计的话那不就不准了，那我问你，你的插值算法准吗？你要准，那你就自己去一点点测量吧。那怎么样根据这些值来估计呢，
　　这就相当于在做线性插值，这里可以放心，因为网格化后的每一个网格的小格在坐标的距离都并不会很大，用什么插值也没有就争议了，这里用一个图来说明
　　

图1　　大家看到了如果说有图那样的情况，是不是那一个等值点就找到了呢，提示大家一下，大家有没有看到这张图里的那个黑点是不是比较靠近10，这可就是线性插值的真实情况了。在一个这样的小格中，这样的点只会有0，2 ，4 个，为什么不会有1，3 个呢，大家想想，如果只有一个，那等值线进去了要怎么样出去呢，3个同样如此，5 个以上那就更不可能，为什么，要注意我们做的是线性插值。这里大家可能想那如果上图中的5变成10，那可怎么办呢，你这样想说明你有点上路了，这个我们会在程序中将其中的一个10加上一个修正值来解决。如果出现有四个等值点的情况，则必须做一下处理，我们就按上图来说明，如果是有四个点，就说明上，右，下边都有等值点，那怎么样就认为8这个点是从那一边出去呢，
　　
　　图2 　　　
　　（3）这种情况是不可能出现的，如果出现，那下一次追踪到这一个小格时就不得不将上下两个点连接起来，这样在画等值线就出现了交叉的情况，这是不允许的，现在要解决的就是究竟是（1）还是（2）的问题，这个问题上也有几种方法，也有比较麻烦的，我只是按照别人的方法用了一个比较简单的方法，就是上下两个点谁8近的方法，如果两个都一样近就看8那面的那个点是偏离上下边哪 个近离上边近就是（1），否则就是（2），按这种方法，上图中就应该是（2）这种方式。
　　要说算法，就有一些不得不说的，那就是算法里用到的算定义的数据结构，
　　Code
struct IsoPoint
　　　　{
　　　　　　public int _column;　　　　 //这里的行列的最大值要比网格的行列数分别小1
　　　　　　public int _row;　　　　　　//因为这里的行列是按行列线（注意是线）来算的
　　
　　　　　　public bool _isHorizon;　　 //等值点是否在X轴上（水平线上）　 true--X　 false--Y
　　　　}

Code　　　　//每条边上的信息（有无等值点，有时点在何处（_rate）
　　　　struct EdgeIsoInfo　　　　
　　　　{
　　　　　　public float _rate;　　 //比率是代表在网格边上的比率位置
　　　　　　public bool _isIsoPoint;　 //在此边上是否有等值点
　　　　}
　　这两种结构就可以将整个网格的各个小格边上是否有等值点，在哪一条边，在哪什么坐标位置（你注意到了_rate了吗？）。
　　在开始前，我们先对整个网格做一次预处理，就是用二维数组保存各个边上的信息 　　　
Code
　　　　privateEdgeIsoInfo[,]_xSide;
　　　　privateEdgeIsoInfo[,]_ySide;
　　其中_xSide, _ySide分别指的是x,y边上的等值点信息。
　　预处理就相当于找到此高程值的所有等值点（为什么，自己想），我这里做的是循环，就是说每一个高程值做一次预处理循环，　有了此高程值相对应的点信息，现在的问题就是如何将所有的这些点分类，分成一条条的等值线，这样就必须要做一些处理，　　那就决定从哪里开始找起，这样就有网格的左，上，右，底边四个边界和网格内部（用于追踪封闭等值线），我们就先做开等值线的追踪，后做封闭等值线的追踪，对于开等值线，我们拿一上面的第一个图来说，也就是追踪左边界，我们从左边界的最下端开始追起，也就是追_ySide[0,0], 看它的那个_isIsoPoints是否为true, 不是就在找_ySide[1,0]，依次类推，假如我们到达10，5 那个位置，也就产图1的情况就达到了示例程序中的从左到右追踪的条件TracingFromLeft2Right（具体见代码），处理完了以后，为了避免再一次
　　追踪到这一点，那我们就将 _isIsoPoints设为false, 这样下次就不会再找到这个点了。
　　对于封闭等值线的处理同样如此，就不再说了。
　　其余详细的看示例代码吧，其实我也是第一次写博客，才知道写东西真的不是那么的容易，以前 看别人一天博客更新一次还嫌慢了，
　　现在到了自己的头上才发现真的不容易，我再说一次，我这里写的，还远远没有将算法讲清楚，也只能想当于介绍了一下子而已，想要
　　弄清楚的大牛们就自己去看源码吧，也可以欢迎和我联系。
　　下一篇就讲光滑了，下次力争写的好一点。

这一次我要说的就是我前面想好要说的等值线的光滑，说到等值线的光滑，我就提一下大家，不知道大家还记不记得在类库中有这样的两个函数，Graphics.DrawCurve 方法和Graphics.DrawBezier方法，说到这两函数，那可能就有人要说，光滑那还不简单，把点传进去不就好了，不错，这样也是一种方法，也达到了比较好的光滑有效果，我在做这个光滑的时候，是先找到的我下面要讲的那种方法，后来才找到这两个函数，那时我可是大笑，可是笑过之后，我还是不得不在这两种方法中做出选择，因为我们要选一个更加的适合我们项目的文法，于是也就有了我下面的这篇文章，至于为什么，我下面在说。
　　等值线的光滑，可能有的比较的了解，她分为两类拟合曲线和逼近曲线，也就是过点型和不过点型，过点型有三次样条和抛物线样条曲线，不过点的有B样条曲线，还有Bezier曲线，当然还有好多，我也不是非常的了解，也象我前面说的，只是做到了够我当前项目用的程度，要了解更深的就自己多查一下这方面的资料，还蛮多的，算法这个东西，我也觉得有时候蛮讨厌的，常常都弄得人头疼，但做出来了又让人很高兴，我做项目最大的感受就是，一个项目比较难做的就是构架和算法，而我在项目里面做的比较多的就是算法，弄得人头疼，当然我说的比较难做也只是说的是开发，至于一个项目里面的什么销售呀什么的，我就没有什么了解了。
　　一种算法的形成都有都其存在的理由，也并不是说这种就一定好，那种就一定不好，就象我用排序算法的时候就一直只用插入排序一样，简单，速度也可以达到要求，我现在就说一种比较简单的过点型（拟合型）的光滑算法，抛物线样条算法，为什么不用三次样条算法，因为三次样条算法有一个要求，就是每一个点都要有导数，这个可是一个比较苛刻的要求，大家要知道，在光滑的时候如果是在B点这样的位置用三次样条算法就不好做了，也可能有人会说那我把X,Y轴对换不就好了吗？

　　那我问大家你把轴对换就一定好了吗，你不信，就多想想，当然也不是说不可以，如果你硬是要那样做的话，你就不得不时时去担心，我的这个点这里要不要将X换成Y，不过我有更好的办法。
　　抛物线样条就是我要说的一类好的方法，这里我也要感谢clever101，因为这个方法就是他做出来的，我只不过是拿他的讲讲而已，抛物线拟合说到底就是三个三个点的拟合，具体的原文大家可以看一下这个：
http://blog.csdn.net/clever101/archive/2006/06/03/771160.aspx
　　我也不多说了，但我还是将他copy到这里 ：
　　＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
　　假如我们采用矢量表达式来表示参数化的二次曲线，那么可以把抛物线的表达式写成如下的一般形式：
　　P(t)=A1+ A2t+ A3t2　　　(0=<t<=1)
　　该抛物线过P1, P2, P3三个点，并且：
　　1.　　　 抛物线以P1点为始点。当参变量t=0时，曲线过P1点；
　　2.　　　 抛物线以P3点为终点。当参变量t=0时，曲线过P3点；
　　3.　　　 当参变量t=0.5时，曲线过P2点,且切矢量等于P3—P1。
　　t=0: P(0)= A1= P1
　　t=1: P(1)= A1 + A2+ A3=P3
　　t=0.5:P(0.5)= A1 + 0.5A2+0.25 A3=P2
　　通过解联立方程，得到三个参数A1 、 A2、 A3分别为：
　　A1 = P1
　　A2=4 P2—P3—3P1
　　A3=2P1+2P3—4P2
　　把求出的这三个系数的值，代入抛物线的表达式P(t)=A1+ A2t+ A3t2

得：P(t)=（2t—3t+1）P1 +(4t—4t2)P2+(4t2—t)P3　 (0=<t<=1)
　　设有一离散型值点列Pi（i=1,2,……,n）,每经过相邻三点作一段抛物线，由于有n个型值点，所以可以做n-2条抛物线段。
　　在这n—2条抛物线段中，第i条抛物线段为经过Pi, Pi+1, Pi+2三点，所以它的表达式应为：Si(ti)=（2t2i—3ti+1）Pi +(4 ti—4 t2i) Pi+1 +(2t2i—ti) Pi+2　　(0=< ti <=1)
　　同理，第i+1条抛物线段为经过Pi+1, Pi+2,Pi+3三点，所以它的表达式应为：Si+1(ti+1)=（2t2i+1—3ti+1+1）Pi+1 +(4 ti+1—4 t2i+1) Pi+2+(2t2i+1—ti+1) Pi+3　　(0=< ti+1 <=1)
　　一般来说，每两段曲线之间的搭接区间，两条抛物线是不可能重合的。如下图所示：

　　显然，对于拟合曲线来说，整个型值点必须只能用一条光滑的曲线连接起来。为了做到这一点，必须找一种方法把Si和Si+1 这样的曲线段的共同区间结合起来。这种方法就是加权合成方法。
　　我们设共同区间的函数是Pi+1(t)=f (T ) Si(ti)+g ( T) Si+1(ti+1). 其中f (T ) 和 g ( T) 是权函数。在抛物样条曲线中我们取简单的一次函数为权函数，且具有互补性，设
　　f (T ) =1—T
　　g ( T) =T
　　这样Pi+1(t)= (1—T ) Si(ti)+ T Si+1(ti+1).因为 函数中有T、ti和ti+1三个参数，因此接下来我们的工作是统一参数。
　　我们可以三个参变量统一形式为：
　　T=2t
　　ti=0.5+t
　　ti+1=t
　　这样
　　Pi+1(t)= （—2t3+4t2—t）Pi +(12t3—410t2+1) Pi+1 +(—12t3+8t2+t) Pi+2 +(4t3—2t2) Pi+3 　　(0=< ti <=0.5)
　　从几何意义上说，函数Pi+1(t)表示的上图的点Pi+1,到Pi+2 之间的线段。但是我们应该看到这种方法从n个点中只能得到n—3段曲线。但是n个型值点应有n—1段曲线。一个直接的想法是添加两个辅助点。那么如何添加呢？
　　方法一：两个辅助点为P0和Pn+1,P0=P1，Pn+1= Pn ,这样画出的曲线为一条不闭合的自由曲线。
　　方法二：添加三个辅助点，P0、Pn+1和Pn+2，然后P0=Pn，Pn+1= P1, ，Pn+2= P2,这样画出的曲线为一条闭合的曲线。
　　＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
　　这种方法比较好的解决了我的问题，这种方法一，它过点，过所有的观察点，二，也没有了什么X,Y轴对换的问题，三，相对简单，不象别的程序算法都是一大堆一大堆的，现在我也没有发现它有什么问题，这也就回答了我为什么用这种方法了。顺便说一句，这个也只是一个思路而已，大家要自己多想想，看有没有不适合自己的地方，自己修改吧，我开始的时候完完全全用的是他的可是没有达到我的要求，我也是做了一些小改动的。
　　附图两张：

　　附上示例程序：以后我会将我做的一整套等值线都上传的，这里也感谢clever101。

现在继续我的等值线追踪算法系列中的最后一个话题，等值线的填充。
　　等值线一系列的问题中，严格上讲，真正是我做的还只有我要讲到的填充了，我前面也说过，算法是比较会让人头疼的一种东西，等值线的填充也是，前面的那两个追踪和光滑多多少少都有文章啦，源码呀来给我参考，虽然说并没有做的比较好，但还是有个参考做的还基本达到了项目的要求了，但这个填充呢，资料非常的少，更不要说源码了，我做出来的这个例子其实只是一个示例，它还有一些问题没有解决，就是颜色的梯度变化，比如说中间填了60这个颜色，上面究竟填55，还是填65的问题，就没有解决，当时做完了这个示例，老师就让我去做三维的控件了，也就没有深入的解决这个问题，如果网友谁做好了，不知道可否发我一分，大家交流交流。
　　说到三维控件，我就多说几句，三维这方面的资料也象那些什么样的算法一样少的可怜，最烦人的就是没有人教，学起来非常的痛苦，记得其中有两个问题，一个就是平移的问题，就是三维空间中的物体要鼠标移动到哪里，物体就移动到哪里，这个问题在我学习MDX的两个月中一直都没有做好，这当中我也问了好多人，别人没有一个人回答的，记得只有CSDN中有一个网友提供了一个方法，思路比较好，可是还是没有解决问题，这个问题还是后来我无意中看到浙大老师的一个示例时才想到的，没有想到非常的简单，只要一句代码就好了。所以说有人教的话，也不会拖两个多月了，还有一个问题就是旋转，也是一个烦人的问题，这个可不是说一句代码就可以做好的，也是做了好久才做好的，以后我也写一个什么系列，只写这两个问题，就平移和旋转，因为我们这学Directx时这两个问题烦了人两个多月，我一直想在别人游戏公司这两个问题肯定不是问题，可是别人就是不告诉你，唉。
　　回到主题上来，做等值线方面的人一定会有“等值线生成与填充算法”，孙桂茹写的那篇论文，我的填充算法就是根据那个做出来的，当然我也找了一些这方面资料对比以后才用这个方法的，不过也好，这样我的查阅资料的能力也提高了不少，越是资料少就越对人的查阅能力要求更高，没有人教就越对人的自学能力要求更高。在那篇论文中，她的算法思路大概就是：

　　等值线只会有上面的a), b), c), d)这四类，它们的覆盖关系为:
　　a 区内部决不会出现b 区,
　　b 区内部决不会出现c 区。
　　因此, 在填充时只要依据一定的顺序依次填充第三(c)、第二(b)、第一种等值线(a)与网格边界所围的区域， 然后按照由外层向内层的顺序填充第四种等值线(d)所围的区域就可以完成对整个区域的填充。
　　整个算法的基本描述如下:
　　1) 按起点纵坐标从下至上的顺序对起点在左边界上的等值线排序;
　　2) 按起点横坐标从左至右的顺序对起点在上边界上的等值线排序;
　　3) 按起点纵坐标从上至下的顺序对起点在右边界上的等值线排序;
　　4) 按起点横坐标从右至左的顺序对起点在下边界上的等值线排序;
　　5) 按起点横坐标从左至右的顺序对内部封闭的等值线排序;
　　6) 填充第三种等值线与网格下边界或左边界以及起点和终点所在的边界所围的区域。对最后一条等值线, 则还需填充与网格上边界或右边界所围的区域;
　　7) 填充第二种等值线与起点和终点所在的边界以及这二边界相交的顶点所围的区域;
　　8) 填充第一种等值线与起点和终点所在的边界的顶点所围的区域;
　　9) 填充内部封闭等值线所围的区域.
　　按着这个做下去就可以达到目的，算法无论别人怎么讲，你都是学不会的，只有你自己悟到了，才能算自己真正的学会了，我这里就不贴代码了。但我告诉大家一个比较好的方法，就是找一个要填充的等值图，按照上面给的方法，在纸上把它们里面的规律找出来，然后自己写程序就可以了。这个我现在没有讲清楚是因为过两天老师回来以后，我就会去更加深入的做一下，前面也说了，我只做了一个示例，后来我就没有做这个了，而由老师接手，想把这个集成到我们的项目里面去，可是时间太紧，别人不停地催老师小结报告，老师没有办法，就只好将这个功能推迟，现在老师回来了，也就到了时候将这个深入下去了，当我做好了以后，我会把整套的程序完整的上传上来，也将这篇文章重新写一下。