等值线的追踪算法

转载自：http://dev.21tx.com/2008/12/15/13505.html

        在等值线的追踪中有一类相对来说比较简单的，也就是涉及不是很深的问题，那就是本文要说的基于规则的等值线的追踪网格。先说一下等值线的形成，等值线的绘制是有两种方法的，一种就是完全测量，也就是在实际的施工中找到所有的等值点，标记上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）这种方式。要说算法，就有一些不得不说的，那就是算法里用到的算定义的数据结构，

struct IsoPoint
{
　 public int _column;　　　　 //这里的行列的最大值要比网格的行列数分别小1
　 public int _row;　　　　　　//因为这里的行列是按行列线（注意是线）来算的
　 public bool _isHorizon;　　 //等值点是否在X轴上（水平线上）　 true--X　 false--Y
}

//每条边上的信息（有无等值点，有时点在何处（_rate）
struct EdgeIsoInfo　　　　
{
　　public float _rate;　　 //比率是代表在网格边上的比率位置
　　public bool _isIsoPoint;　 //在此边上是否有等值点
}

这两种结构就可以将整个网格的各个小格边上是否有等值点，在哪一条边，在哪什么坐标位置（你注意到了_rate了吗？）。

　　在开始前，我们先对整个网格做一次预处理，就是用二维数组保存各个边上的信息 　　　

privateEdgeIsoInfo[,]_xSide;
privateEdgeIsoInfo[,]_ySide;

其中_xSide, _ySide分别指的是x,y边上的等值点信息。预处理就相当于找到此高程值的所有等值点（为什么，自己想），我这里做的是循环，就是说每一个高程值做一次预处理循环，有了此高程值相对应的点信息，现在的问题就是如何将所有的这些点分类，分成一条条的等值线，这样就必须要做一些处理，那就决定从哪里开始找起，这样就有网格的左，上，右，底边四个边界和网格内部（用于追踪封闭等值线），我们就先做开等值线的追踪，后做封闭等值线的追踪，对于开等值线，我们拿一上面的第一个图来说，也就是追踪左边界，我们从左边界的最下端开始追起，也就是追_ySide[0,0], 看它的那个_isIsoPoints是否为true, 不是就在找_ySide[1,0]，依次类推，假如我们到达10，5 那个位置，也就产图1的情况就达到了示例程序中的从左到右追踪的条件TracingFromLeft2Right（具体见代码），处理完了以后，为了避免再一次追踪到这一点，那我们就将 _isIsoPoints设为false, 这样下次就不会再找到这个点了。对于封闭等值线的处理同样如此，就不再说了。

　　其余详细的看示例代码吧，其实我也是第一次写博客，才知道写东西真的不是那么的容易，以前 看别人一天博客更新一次还嫌慢了，现在到了自己的头上才发现真的不容易，我再说一次，我这里写的，还远远没有将算法讲清楚，也只能想当于介绍了一下子而已，想要弄清楚的大牛们就自己去看源码吧，也可以欢迎和我联系。下一篇就讲光滑了，下次力争写的好一点。