posts - 116,  comments - 267,  trackbacks - 1

转载请注明作者及出处,谢谢

这两天手头有个项目,需要绘制等值线,本以为是一个很简单的事情,没有想到刚开始就发现竟然无从着手,研究了一个星期,终于把线条画出来了,基本思路是先三角网剖分,然后再等值线追踪,最后绘制;没有对等值线进行光滑处理,示例图中看起来比较光滑是因为取点比较密集,也没有打算进行等值线填色,因为项目中没有这个需求,(而且在我的项目中高程点是网格状分布,而不是离散点,因此我做的三角网剖分简单,但是等值线追踪算法是完全满足离散点要求的)。先上几个效果图:

示例图(黄颜色圆圈代表光源,高程值为光源照度)

图1

图2

图3

等值线标注示意图

效果一:高程值压线了

效果二:高程值在线条下方

效果三:高程值没有按照等值线方向旋转

效果四:在很小的封闭式等值线圈内进行标注

绘制等值线的理论有很多,但是通常大家采用的是三角网剖+等值线追踪,网上有很多论文可以参考:

Delaunay三角网剖分算法

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个点取一点。

有了三角网,接下来应该去进行等值线追踪了。

稍事休息,接下来讲等值线追踪。。。

posted on 2011-04-29 14:33  think8848  阅读(...)  评论(...编辑  收藏