[原]C#繪製等值線一 基本概念及三角網剖分

轉載請註明作者及出處，謝謝

這兩天手頭有個項目，需要繪製等值線，本以為是一個很簡單的事情，沒有想到剛開始就發現竟然無從著手，研究了一個星期，終於把線條畫出來了，基本思路是先三角網剖分，然後再等值線追蹤，最後繪製；沒有對等值線進行光滑處理，樣本圖中看起來比較光滑是因為取點比較密集，也沒有打算進行等值線填色，因為項目中沒有這個需求，(而且在我的項目中高程點是網格狀分布，而不是離散點，因此我做的三角網剖分簡單，但是等值線追蹤演算法是完全滿足離散點要求的)。先上幾個：

樣本圖(黃顏色圓圈代表光源，高程值為光源照度)

圖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個點取一點。

有了三角網，接下來應該去進行等值線追蹤了。

稍事休息，接下來講等值線追蹤。。。