转载请注明作者及出处，谢谢
这两天手头有个项目，需要绘制等值线，本以为是一个很简单的事情，没有想到刚开始就发现竟然无从着手，研究了一个星期，终于把线条画出来了，基本思路是先三角网剖分，然后再等值线追踪，最后绘制；没有对等值线进行光滑处理，示例图中看起来比较光滑是因为取点比较密集，也没有打算进行等值线填色，因为项目中没有这个需求，(而且在我的项目中高程点是网格状分布，而不是离散点，因此我做的三角网剖分简单，但是等值线追踪算法是完全满足离散点要求的)。先上几个效果图：
示例图(黄颜色圆圈代表光源，高程值为光源照度)
图1
图2
图3
等值线标注示意图
效果一：高程值压线了
效果二：高程值在线条下方
效果三：高程值没有按照等值线方向旋转
效果四：在很小的封闭式等值线圈内进行标注
绘制等值线的理论有很多，但是通常大家采用的是三角网剖+等值线追踪，网上有很多论文可以参考：
Delaunay三角剖分(Delaunay Triangulation)相关知识
这里是一个绘制、填充等值线的源代码(《C#实现基于三角网的等值线追踪及填充算法》）
这里是一些三角网剖分源代码，有各种语言的实现
这里是一些其他等值线绘制方法
这里是论文的链接地址，该链接包含很多论文地址。
基本上来说，在《三角网格等值线自动生成方法及程序实现》一文中，提供的理论就足以指导绘制作业了，但是不幸的是，相当多的程序员可能不会一下子就能行成一个比较清晰明确的思路，我就是其中之一，看了好多的文章才明白，噢！原来是这么回事啊。另外还是那句老话：“眼里过千遍，不如手里过一遍”，所谓纸上得来终觉浅，绝知此事要躬行，看算法，看代码，可能比较难以一下子明白其中的原理，很多时候人家用了一个算法(或数据结构)，你可能很不以为然，但是当你自已动手时，才发现原来那样用有那样用的妙处(或只能那么用)，眼高手低永远是我等程序员的通病。
ok，进入正题
在我的项目中，是计算在一个空间内部，假设排布了n盏灯，那么需要计算出空间内部每一个点上的照度是多少，同时绘制出照度的等值线来，计算照度不在本文计论范围之内，于是我们假设把空间划分成了一个矩阵，矩阵中每个点上的照度值均已获得。于是很自然想到“使用网格绘制等值线”这么一个关键词，发现居然很少有相关资料，于是只能把把关键词修改为“绘制等值线”，于是发现了很多关于绘制等值线方法的论文，但是基本上前面都带了另一个词：三角网。为什么要三角网呢？我们先来看一下使用网格绘制时一种情形
图1
在使用网格绘制等值线时，假定我们是从这个网格的左边起开始查找等值线的走向，则等值线有三个方向可以“出去”，判断和确定等值线的走向将是一个非常困难的事情，而且极有可能出现锯齿状等值线(如下图所示)
图2
而使用三角网为基础来绘制等值线时，则是下图的情况
图3
在这里我有一个问题，到目前为止还没有解决，但是也没有找到解决问题的理论依据，在使用网格时，需要判断等值线的走向，有个论文中提出一个原则：假设等值线的“入口”在左边上(图1A点)，则在网格中，如果等值线的“出口”有多个(图1B,C,D点)，那么首先要看到等值线到A的趋势更倾向于择哪个点(B,C,D)，其次是看哪个点到A点的距离近则取哪个。但是实际上这样做的代价太大了，任何人要去写一段判断趋势的代码，肯定都是个头疼的问题，因此实际上的作法是：在由左向右的追踪过程中，先看哪个点的X值和A点的X值接近则取哪一个，如果有多个点的X值和A点的X值相等，则看哪个点到A点的距离小则取哪一个，如果这个距离也相等，则只能通过趋势判断了(见《基于规则网格等值线生成算法及其应用》1.2.2)。实际从程序员的角度来讲，只要你明确表示在程序中会用到某个功能，那么你一年用n+1次和用1次是没有区别的，因此如果你想把程序写的能覆盖所有的可能，判断趋势看起来是不可避免的。
虽然三角形简化了上述这种在网格中追踪等值线的难度，但是并没有完全解决需要判断趋势的问题，为啥咧？使用三角形时只可能有两个“出口”，使用网格中的原则基本上就可以判断取哪个点了，但是理论上在点A所在的边上，仍一个点，到B，和C的X值和距离是相等的，见下图
图4
在图4中，三角形ABC是一个等腰三角形，这就意味着A点到B以及C的距离相等，而且B和C的X值也相等...理论上来讲，如果A是起点，哪还谈什么趋势呢？还要不要让人活了...
OK，现在让我们一起默念“阿门”，求助于概率吧，不要出现这种情况，反正我是不打算处理这种情况的，出了Bug再说吧。值的说明的是，在约大多数的论文中，没有提到该如何去处理这情况，有个别论文中指出：“每个三角形上不可能三边都有同值的等值点，另一边上必定有同值的等值点”(见《如何根据离散点自动绘制等值线(等高线)之 三角形法》1.b)，但是这个说法我持怀疑态度，但我现在不打算处理这种情况(实际上，我在代码中连距离和X值都没有处理)。
有了理论作为依据，接下来的事情就是先把矩阵转换为三角网，怎么转换是个问题：实际上很多时候关于Delaunay的论文都是在假定数据(如灯具照度、高度、降水量等)是离散数据，即不是规则分布的，因此在这些论文中都是在讲如果把这些离散点划分为三角形，当然，能把离散点剖分出三角网，则么规则点就不在话下了。但在我的项目中是不存在这些需求，因此我需要做的就是把矩阵中的每个网格一分为二，切成两个直角三角形，一率从左上角向右下角切，如图所示：
图5
切的有点密，这个可以视具体情况自已定夺。
在这里，需要把三角形、边、点的数据结构给出来：
VPoint代表一个点，由X，Y轴坐标及高程值构成
public class VPoint
{
public float X { get; set; }
public float Y { get; set; }
public decimal Value { get; set; }
public override bool Equals(object obj)
{
if (obj is VPoint)
{
var tmp = obj as VPoint;
return this.X == tmp.X && this.Y == tmp.Y;
}
return false;
}
public override int GetHashCode()
{
return base.GetHashCode();
}
}
Edge代表一个边，由两个点构成，这里需要稍作解释，在图5中，有一些三角形的一条边是整个三角网的边框，即空间的四个边，开放式的等值的起点和终点必然会落在这些最外层的边上，因此需要对这些最外层边作一个标记，如何标记呢？很显然，三角网中除了这些最外层边，其他所有的边都会同时属于两个三角形，因此，我们的Edge对象有两个属性：Trangle1和Trangle2，这两个属性分别标识这个Edge对象所属的两个三角形，那么这些最外层的Edge对象比较可怜，他们只属于一个Trangle,属性Trangle2的值一定是null，参见Trangle类型。
public class Edge
{
public VPoint P1 { get; set; }
public VPoint P2 { get; set; }
public Triangle Triangle1 { get; set; }
public Triangle Triangle2 { get; set; }
public bool IsBorder
{
get
{
return Triangle2 == null;
}
}
public override bool Equals(object obj)
{
if (obj is Edge)
{
Edge tmp = obj as Edge;
return this.P1 == tmp.P1 && this.P2 == tmp.P2;
}
return false;
}
public override int GetHashCode()
{
return base.GetHashCode();
}
}
Trangle代表一个三角形，由三个边构成，我总是把水平的那条边称为A，竖直的那条边称为B，斜边称为C，其他的属性值以后会用，这里就不一一介绍了，当然，现在代码还没有进行最终整理，有可能有些属性也许会被干掉，最终只留下必要的。
public class Triangle
{
private Edge a, b, c;
/// <summary>
/// 邻边
/// </summary>
public Edge A
{
get
{
return this.a;
}
set
{
this.a = value;
if (this.a.Triangle1 == null)
{
this.a.Triangle1 = this;
}
else
{
this.a.Triangle2 = this;
}
}
}
/// <summary>
/// 对边
/// </summary>
public Edge B
{
get
{
return this.b;
}
set
{
this.b = value;
if (this.b.Triangle1 == null)
{
this.b.Triangle1 = this;
}
else
{
this.b.Triangle2 = this;
}
}
}
/// <summary>
/// 斜边
/// </summary>
public Edge C
{
get
{
return this.c;
}
set
{
this.c = value;
if (this.c.Triangle1 == null)
{
this.c.Triangle1 = this;
}
else
{
this.c.Triangle2 = this;
}
}
}
/// <summary>
/// 是否使用,-1：临时使用; 0:未使用; 1:使用
/// </summary>
public int UsedStatus { get; set; }
public Edge BorderEdge
{
get
{
if (this.A.IsBorder) { return this.A; }
if (this.B.IsBorder) { return this.B; }
return null;
}
}
public Edge[] Get2OtherEdge(Edge relativeEdge)
{
Edge[] edges = new Edge[2];
if (relativeEdge == this.A)
{
edges[0] = this.B;
edges[1] = this.C;
}
else if (relativeEdge == this.B)
{
edges[0] = this.A;
edges[1] = this.C;
}
else
{
edges[0] = this.A;
edges[1] = this.B;
}
return edges;
}
public VPoint FindPoint(float x, float y)
{
if (this.A != null)
{
if (this.A.P1.X == x && this.A.P1.Y == y) { return this.A.P1; }
if (this.A.P2.X == x && this.A.P2.Y == y) { return this.A.P2; }
}
if (this.B != null)
{
if (this.B.P1.X == x && this.B.P1.Y == y) { return this.B.P1; }
if (this.B.P2.X == x && this.B.P2.Y == y) { return this.B.P2; }
}
if (this.C != null)
{
if (this.C.P1.X == x && this.C.P1.Y == y) { return this.C.P1; }
if (this.C.P2.X == x && this.C.P2.Y == y) { return this.C.P2; }
}
return null;
}
public Edge FindEdge(VPoint p1, VPoint p2)
{
if (this.A != null)
{
if (this.A.P1 == p1 && this.A.P2 == p2) { return this.A; }
}
if (this.B != null)
{
if (this.B.P1 == p1 && this.B.P2 == p2) { return this.B; }
}
if (this.C != null)
{
if (this.C.P1 == p1 && this.C.P2 == p2) { return this.C; }
}
return null;
}
}
OK，现在是时候说说我是如何建立一个IList<Trangle>的吧，即生成三角网列表
一开始，需要解决“有木有”的问题，于是我先建立一个IList<Trangle>的实例result，然后生成四个VPoint对象，代表一个网格的四个角，然后再生成五个Edge对象，再使用五个Edge生成两个Trangle对象，然后把两个Trangle对象加入到result中，以后每次都去查找result中有没有指定X,Y的VPoint，以及查找result中有没有指定P1和P2的Edge，如果有的话，就分别拿出来去构成新的Edge（对于VPoint来说）或者Trangle（对于Edge来说），如果没有的话就new一个实例出来，然后再通过Trangle对象存到result中，参见Trangle.FindPoint方法和Trangle.FindEdge方法，因为大量使用Linq方法查询数据，造成的结果是45 × 27的矩阵，生成三角网需要花费00:00:04.7031250！后来我改进了方法，使用一个二维数组来保存VPoint对象实例(VPoint[x, y])，使用一个四维数组来保存Edge对象实例(Edge[p1x, p1y, p2x, p2y])，直接用矩阵的序号来获取可能已存在的实例，果然性能得到极大的提升，运行耗时仅00:00:00.0156250。
for (int x = 0; x < xCount - step; x += step)
{
for (int y = 0; y < yCount - step; y += step)
{
var p1 = matrix[x, y];
var p2 = matrix[x + step, y];
var p3 = matrix[x + step, y + step];
var p4 = matrix[x, y + step];
VPoint tmpP1, tmpP2, tmpP3, tmpP4;
Edge edge1, edge2, edge3, edge4, edge5;
Triangle triangle1 = new Triangle(), triangle2 = new Triangle();
tmpP1 = this.FindOrCreateNewPoint(tmpMatrix, x, y, p1, zoomFactor);
tmpP2 = this.FindOrCreateNewPoint(tmpMatrix, x + 1, y, p2, zoomFactor);
tmpP3 = this.FindOrCreateNewPoint(tmpMatrix, x + 1, y + 1, p3, zoomFactor);
tmpP4 = this.FindOrCreateNewPoint(tmpMatrix, x, y + 1, p4, zoomFactor);
edge1 = this.FindOrCreateNewEdge(tmpEdges, tmpMatrix, x, y + 1, x + 1, y + 1);
edge2 = this.FindOrCreateNewEdge(tmpEdges, tmpMatrix, x, y, x, y + 1);
edge3 = this.FindOrCreateNewEdge(tmpEdges, tmpMatrix, x, y, x + 1, y);
edge4 = this.FindOrCreateNewEdge(tmpEdges, tmpMatrix, x + 1, y, x + 1, y + 1);
edge5 = this.FindOrCreateNewEdge(tmpEdges, tmpMatrix, x, y, x + 1, y + 1);
triangle1.A = edge1;
triangle1.B = edge2;
triangle1.C = edge5;
triangle2.A = edge3;
triangle2.B = edge4;
triangle2.C = edge5;
result.Add(triangle1);
result.Add(triangle2);
}
}
matrix是原始数据，在我的项目是照度点矩阵，step的作用是由原始数据中生成三角网时可以跳过step个点取一点。
有了三角网，接下来应该去进行等值线追踪了。
稍事休息，接下来讲等值线追踪。。。
