等值線的追蹤演算法

來源:互聯網
上載者:User

   在等值線的追蹤中有一類相對來說比較簡單的,也就是涉及不是很深的問題,那就是本文要說的基於規則的等值線的追蹤網格。
  先說一下等值線的形成,等值線的繪製是有兩種方法的,一種就是完全測量,也就是在實際的施工中找到所有的等值點,標記上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) 填充內部封閉等值線所圍的地區.
  按著這個做下去就可以達到目的,演算法無論別人怎麼講,你都是學不會的,只有你自己悟到了,才能算自己真正的學會了,我這裡就不貼代碼了。但我告訴大家一個比較好的方法,就是找一個要填充的等值圖,按照上面給的方法,在紙上把它們裡面的規律找出來,然後自己寫程式就可以了。這個我現在沒有講清楚是因為過兩天老師回來以後,我就會去更加深入的做一下,前面也說了,我只做了一個樣本,後來我就沒有做這個了,而由老師接手,想把這個整合到我們的項目裡面去,可是時間太緊,別人不停地催老師小結報告,老師沒有辦法,就只好將這個功能延遲,現在老師回來了,也就到了時候將這個深入下去了,當我做好了以後,我會把整套的程式完整的上傳上來,也將這篇文章重新寫一下。 

聯繫我們

該頁面正文內容均來源於網絡整理,並不代表阿里雲官方的觀點,該頁面所提到的產品和服務也與阿里云無關,如果該頁面內容對您造成了困擾,歡迎寫郵件給我們,收到郵件我們將在5個工作日內處理。

如果您發現本社區中有涉嫌抄襲的內容,歡迎發送郵件至: info-contact@alibabacloud.com 進行舉報並提供相關證據,工作人員會在 5 個工作天內聯絡您,一經查實,本站將立刻刪除涉嫌侵權內容。

A Free Trial That Lets You Build Big!

Start building with 50+ products and up to 12 months usage for Elastic Compute Service

  • Sales Support

    1 on 1 presale consultation

  • After-Sales Support

    24/7 Technical Support 6 Free Tickets per Quarter Faster Response

  • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.