專利名稱：基于重力場模型的管道狀物體拉直的圖像處理方法
技術領域：
本發明涉及的是一種圖像處理技術領域的方法，具體是一種基于重力場模型 的管道狀物體拉直的圖像處理方法。
背景技術：
虛擬可視化技術是一種觀察管道狀物體的三維圖像可視化技術。這個技術的 優勢在于無需破環或接觸物體本身，通過計算機圖像處理模擬切開、展平、外翻 等方式，以觀察管道內壁結構和狀態。虛擬拉直可以將彎的管道映射到直的柱面 上，有助于虛擬展平、外翻可視化技術的實施，既可以保持轉彎處的細節信息， 又能夠提供直觀的整體信息。
經過對現有技術的文獻檢索發現，美國Iowa大學放射系王革(Ge Wang)教授 等于1995年首次提出了基于電場模型和中心路徑拉直的曲截面結腸拉直展平方 法，闡述了結腸拉直在虛擬可視化技術中的重要作用(G. WangandM. W. Vannier， 〃GI tract unraveling by spiral CT, 〃 in Proc. SPIE， San Diego, CA, USA, 卯.307-315, 1995)。這個方法的不足在于對于管道狀物體彎曲較大(中心路 徑曲率較大)的區域，由于同種電荷的互斥作用，基于電場模型的曲截面呈現不 對稱的分布，這將導致拉直變形后在該區域產生嚴重的失真；對于管道內部由于 凸起、凹陷等原因造成的不對稱結構區域，基于中心路徑拉直的方法會使管道狀 物體產生不正常的扭曲導致失效，從而無法得到反映內腔的客觀準確的三維拉直 圖像。
發明內容
本發明針對現有技術存在的上述不足，提供一種基于重力場模型的管道狀物 體拉直的圖像處理方法，對于管道狀物體彎曲較大的區域，以及管道內部由于凸 起、凹陷等原因造成的不對稱結構區域，能夠進一步獲得優化的管道狀物體拉直 后的圖像，避免了現有技術中基于電場模型和中心路徑拉直方法所造成的拉直算 法失效問題。
本發明是通過以下技術方案實現的，本發明包括以下步驟第一步、通過計算機斷層掃描成像法、工業計算機斷層掃描成像法、核磁共 振成像法或超聲成像法獲得原始三維圖像數據，然后對原始三維圖像數據采用濾 波去噪法、對比度調整法或插值法進行優化預處理，最后通過區域增長或水平集 方法進行圖像分割，得到管道狀物體的內壁三維圖像數據和外壁三維圖像數據。
所述的內壁三維圖像數據是指通過區域增長或水平集方法獲得的管道狀物 體的內部空腔三維圖像數據。
所述的外壁三維圖像數據是指利用獲取的內壁三維圖像數據，通過膨脹腐 蝕或水平集方法獲得的管道狀物體的內部空腔三維圖像數據、管道狀物體的管壁 三維圖像數據的總和。
第二步、將外壁三維圖像數據劃分為結構化網格、掃略網格或自由網格中的 一種并生成外壁三維有限元模型。
所述的劃分是指通過Medial Axis算法或Advancing Front算法進行劃分處理。
所述的有限元網格的網格單元為四面體單元或六面體單元。
第三步、對外壁三維有限元模型施加邊界條件和重力場載荷。
所述的邊界條件是指在外壁三維有限元模型的一端取一橫截面，設定位于該 橫截面上的所有節點的位移和速度均為0。
所述的重力場載荷是指設定垂直于施加邊界約束的橫截面，且指向外壁三維 有限元模型內側的方向為重力場方向。
第四步、依據所施加的邊界約束和重力載荷對外壁三維有限元模型進行非線 性有限元變形計算求解，生成拉直變形后的外壁三維有限元模型。
所述的非線性有限元變形計算求解，是指列出虛功方程，使用Newton-Raphson方法、修正的Newton-Raphson方法、準Newton方法等迭代算法中的一種 計算虛功方程的極值點。
所述的虛功方程是指:變形體中任意滿足平衡的力系在任意滿足協調條件的 變形狀態上作的虛功等于零，即體系外力的虛功與內力的虛功之和等于零。
所述的協調條件是指在變形體內部連續可導。
第五步、使用基于特征點的非剛體配準的方法實現內壁三維圖像數據的拉直 變形將外壁三維有限元模型上的各節點作為源特征點(即內壁三維圖像數據的特征點)，將第四步中計算求解得到的拉直變形后的外壁三維有限元模型上的各節 點作為目標特征點(即拉直變形后的內壁三維圖像數據的特征點)，根據源特征點 與目標特征點的對應關系，應用非剛體配準方法，得到內壁三維圖像數據中任意 一點到拉直變形后的內壁三維圖像數據中對應點的映射關系，從而得到拉直變形 后的內壁三維圖像數據。
所述的非剛體配準方法是指薄盤樣條法、近似薄盤樣條法、體積樣條法或 彈性體樣條法中的一種。
第六步、采用中心路徑生成算法計算外壁三維圖像數據的中心路徑和拉直變 形后的外壁三維圖像數據的中心路徑。
所述的中心路徑生成算法是指根據拓撲學細化法、距離變換法、水平集方法 或手工標定法得到中心路徑的方法。
所述的拉直變形后的外壁三維圖像數據是指利用獲取的拉直變形后的內壁 三維圖像數據，通過膨脹腐蝕或非剛體配準等方法中的一種獲得的拉直變形后的 管道狀物體的內部空腔三維圖像數據、拉直變形后的管道狀物體的管壁三維圖像 數據的總和。
所述的中心路徑是指管道狀物體內部的一根中軸軌跡，該中軸軌跡上各點都 滿足最優中心點判定準則。
第七步、基于第五步中得到的特征點進行旋轉失真校正處理根據第五步中 得到的源特征點與目標特征點的對應關系，計算拉直變形后的內壁三維圖像數據 相鄰數據層的旋轉角度，對拉直變形后的內壁三維圖像數據第二層以后的每一層， 進行旋轉變換，得到旋轉失真校正后的內壁三維圖像數據。
第八步、基于第五步中得到的特征點和第七步的結果，進行長度失真校正處 理，將長度失真校正后的內壁三維圖像數據重新繪制，以三維形式顯示出來。觀 察視角位于管道內部或外部，并可由用戶任意選取。
長度失真校正的方法包括重采樣方法和基于特征點的非剛體配準方法。
所述的重采樣方法是指根據第五步中得到的源特征點與目標特征點的對應 關系，保持外壁三維圖像數據的中心路徑與長度失真校正后的外壁三維圖像數據 的中心路徑長度相同，對旋轉失真校正后的內壁三維圖像數據重采樣，生成長度 失真校正后的內壁三維圖像數據。所述的基于特征點的非剛體配準方法是指根據外壁三維圖像數據的中心路 徑與長度失真校正后的外壁三維圖像數據的中心路徑長度相同的條件，重新計算 第五步中得到的源特征點與目標特征點的對應關系，生成用于長度失真校正的特 征點及對應關系，并據此應用非剛體配準方法，得到旋轉失真校正后的內壁三維 圖像數據中任意一點到長度失真校正后的內壁三維圖像數據中對應點的映射關 系，生成長度失真校正后的內壁三維圖像數據。
本發明相對于現有的管道狀物體拉直的圖像處理方法的有益效果在于1、對 于管道狀物體彎曲較大(中心路徑曲率較大)的區域，避免了由于同種電荷的互 斥作用使得基于電場模型的曲截面呈現不對稱的分布，導致管道狀物體拉直后在 該區域產生嚴重失真。2、對于管道內部由于凸起、凹陷等原因造成的不對稱結構 區域，避免了基于中心路徑拉直的方法會使管道狀物體產生不正常的扭曲導致失 效的問題。
圖1為實施例原始三維圖像數據。
圖2為實施例采用現有技術處理后示意圖。
圖3為實施例采用基于重力場模型的拉直方法變形后示意圖。
具體實施例方式
下面對本發明的實施例作詳細說明，本實施例在以本發明技術方案為前提下 進行實施，給出了詳細的實施方式和具體的操作過程，但本發明的保護范圍不限 于下述的實施例。
本實施例具體通過以下步驟實施
第一步、對己獲得的一個管道狀物體CT斷層圖像數據進行濾波處理，圖1所 示，使用三維區域生長算法從斷層圖像中分割出粗略的管道狀物體內部空腔。再 使用三維水平集分割算法，將區域生長算法的分割結果作為初值，將原斷層圖像 作為速度圖像，進一步分割出精確的物體內部空腔，即內壁三維圖像數據。使用 方向梯度方法從原斷層圖像中提取出管道狀物體的內外表面輪廓，此輪廓作為速 度圖像，內壁三維圖像數據作為初值，再次使用三維水平集分割算法提取出物體 的外壁三維圖像數據。
第二步、對外壁三維圖像數據劃分有限元網格，生成外壁三維有限元模型。選擇四面體網格單元劃分結構化網格。賦值單元的密度屬性、楊氏模量、泊松比。 其中泊松比可設為0，以保證拉直過程不會影響徑向寬度。
第三步、對于外壁三維有限元模型的一端，取一橫截面，位于截面上的所有 節點施加ENCASTRE約束，即定義其位移和速度均為0 。施加全局重力場，重力加 速度取9.8m"2，場方向為垂直于施加邊界約束的截面指向模型內側的方向。
第四步、對外壁三維有限元模型進行非線性有限元變形計算求解。設管道狀 物體的占有域為/2，其邊界為廠，對于參考構型中的材料點X，虛功(率)方程 為<formula>formula see original document page 10</formula>
其中，^ 為總虛功率，^mt為內部虛功率，^^'為由物體外力b(X,0和指定的
面力i(X,O所引起的外部虛功率，p(X,O為密度，"③為維數，v(X,O為材料點
的位置矢量的變化率，即速度，tr(X，/)為Cauchy應力，D(Xj)為變形率。使用
Newton-Raphson迭代算法計算虛功方程的極值點，求得拉直變形后的外壁三維有 限元模型及各節點的位置坐標。
第五步、將外壁三維有限元模型上的各節點作為源特征點(即內壁三維圖像 數據的特征點)，將第四步中計算求解得到的拉直變形后的外壁三維有限元模型上 的各節點作為目標特征點(即拉直變形后的內壁三維圖像數據的特征點)，根據源 特征點與目標特征點的對應關系，求得每一對特征點之間的變換向量。再對內壁 三維圖像數據分段利用薄盤樣條插值，得到內壁三維圖像數據中任意一點到拉直 變形后的內壁三維圖像數據中對應點的映射關系，使用基于特征點的非剛體配準 的方法實現內壁三維圖像數據的拉直，得到拉直變形后的內壁三維圖像數據。
第六步、使用距離變換對外壁三維圖像數據的每一點求出該點到外壁三維圖像 數據邊界點的最小距離。采用Dijkstra最短路徑法生成外壁三維圖像數據的中心 路徑即在外壁三維圖像數據中選取中心路徑的起點和終點，然后以外壁三維圖 像數據中每一點到外壁三維圖像數據邊界點的距離值的倒數作為權值，進行最短 路徑搜索，取所有權值最小的點連接起點和終點得到外壁三維圖像數據的中心路 徑。同樣使用距離變換和Dijkstra最短路徑法求得拉直變形后的外壁三維圖像數據的中心路徑。
第七步、基于第五步中得到的特征點和第六步中得到的中心路徑進行旋轉失 真校正處理，具體步驟包括
1 )對于拉直變形后的外壁三維圖像數據的中心路徑上的每一點 C,(/ = 1，2，L ,n)，在拉直變形后的外壁三維有限元模型中選擇與該點距離最小的一
個節點iV,;
2) 平移拉直變形后的外壁三維圖像數據的中心路徑，對于平移后的中心路徑 上的每一點C》、1,2,L ，")，再次在拉直變形后的外壁三維有限元模型中選擇與該 點距離最小的一個節點f;
3) 在外壁三維有限元模型中，可以分別找到唯一一個與W,和《對應的節點
4) 用四組節點中序號相同的節點可以構建兩組向量，5 ^和
UU1 ULI LU UUl
分別計算t/,與下一相鄰向量c/,+i之間的角度《、c/,'與下一相鄰向量
UULt
《+1之間的角度《；
5) 計算角度差
=《-《，z. = l,2,L ,w-1, (2) 作為中心路徑點C,+i所在數據層的旋轉失真的估計；
6) 對拉直變形后的內壁三維圖像數據第二層以后的每一層，進行旋轉變換， 得到旋轉失真校正后的內壁三維圖像數據對拉直變形后的內壁三維圖像數據第 f(/ = l，2,L ，n)層的任意一點(x,力，與之對應的旋轉失真校正后的內壁三維圖像數
據點的坐標(rx,oO滿足
fx, L(x — cx) cos(Zl《一)_ (少一 c_y) sin(zlD + cx」
b，O7 = S i i LL(x - cx) sin(zl《一,)一 (y - c力cos(WM ) + c_y」
其中(CJC,C力為拉直變形后的外壁三維圖像數據的中心路徑在第/層的坐標。
第八步、基于第五歩中得到的特征點、第六步中得到的中心路徑和第七步的
if / = 1， if / > 1，
if,- = 1,
if / > 1，
(3)
(4)
11結果，進行長度失真校正處理，具體步驟包括
1) 利用第七步得到的節點集合iV—^，iV2，L ,iVj計算其中每一點
A^^ = 1,2,L ，w-l)與下一相鄰節點iV^之間的直線段在拉直變形后的外壁三維圖 像數據的中心路徑上的投影長度《；
2) 利用第七步得到的節點集合V^7Vi,iV2，L ,《}計算其中每一點
《(A: = 1，2,L ，n-l)與下一相鄰節點iVw之間的直線段在外壁三維圖像數據的中心 路徑上的投影長度《；
3) 計算第l)步和第2)步所得的投影長度的差
<formula>formula see original document page 12</formula>作為中心路徑點C^所在數據層的長度失真的估計，使用基于特征點的非剛體配
準方法進行長度失真校正對于長度失真校正后的內壁三維圖像數據的第 yt(yt = l，2，L ")層，取點(x,j，A:)作為目標特征點，取點(x，y,m)作為源特征點，其<formula>formula see original document page 12</formula>
其中X和y為拉直變形后的外壁三維圖像數據的中心路徑在第A層的坐標 值。根據源特征點與目標特征點的對應關系，求得每一對特征點之間的變換向量。 再對旋轉失真校正后的內壁三維圖像數據利用薄盤樣條插值，得到旋轉失真校正 后的內壁三維圖像數據中任意一點到長度失真校正后的內壁三維圖像數據中對應 點的映射關系，使用基于特征點的非剛體配準的方法實現長度失真校正，得到長 度失真校正后的內壁三維圖像數據。
第九步、使用VTK等三維可視化工具包將長度失真校正后的內壁三維圖像數 據重新繪制，以用戶可控的方式顯示出來，如圖3所示，用戶可以任意改換視角， 對管道狀物體內部和外部進行觀察。
與現有技術處理方式得到的圖2相比，通過以上步驟，在實現管道狀物體拉 直的同時，解決了高曲率區域拉直后的較大失真和不對稱結構區域導致基于中心 路徑拉直方法失效的問題，能夠方便進行進一步可視化及其他處理操作。
權利要求
1、一種基于重力場模型的管道狀物體拉直的圖像處理方法，其特征在于，包括以下步驟第一步、通過計算機斷層掃描成像法、工業計算機斷層掃描成像法、核磁共振成像法或超聲成像法獲得原始三維圖像數據，然后對原始三維圖像數據采用濾波去噪法、對比度調整法或插值法進行優化預處理，最后通過區域增長或水平集方法進行圖像分割，得到管道狀物體的內壁三維圖像數據和外壁三維圖像數據；第二步、將外壁三維圖像數據劃分為結構化網格、掃略網格或自由網格中的一種并生成外壁三維有限元模型；第三步、對外壁三維有限元模型施加邊界條件和重力場載荷；第四步、依據所施加的邊界約束和重力載荷對外壁三維有限元模型進行非線性有限元變形計算求解，生成拉直變形后的外壁三維有限元模型；第五步、使用基于特征點的非剛體配準的方法實現內壁三維圖像數據的拉直變形將外壁三維有限元模型上的各節點作為源特征點，將第四步中計算求解得到的拉直變形后的外壁三維有限元模型上的各節點作為目標特征點，根據源特征點與目標特征點的對應關系，應用非剛體配準方法，得到內壁三維圖像數據中任意一點到拉直變形后的內壁三維圖像數據中對應點的映射關系，從而得到拉直變形后的內壁三維圖像數據；第六步、采用中心路徑生成算法計算外壁三維圖像數據的中心路徑和拉直變形后的外壁三維圖像數據的中心路徑；第七步、基于第五步中得到的特征點進行旋轉失真校正處理根據第五步中得到的源特征點與目標特征點的對應關系，計算拉直變形后的內壁三維圖像數據相鄰數據層的旋轉角度，對拉直變形后的內壁三維圖像數據第二層以后的每一層，進行旋轉變換，得到旋轉失真校正后的內壁三維圖像數據；第八步、基于第五步中得到的特征點和第七步的結果，進行長度失真校正處理，將長度失真校正后的內壁三維圖像數據重新繪制，以三維形式顯示出來，觀察視角位于管道內部或外部，并可由用戶任意選取。
2、 根據權利要求1所述的基于重力場模型的管道狀物體拉直的圖像處理方法，其特征是，所述的有限元網格的網格單元為四面體單元或六面體單元。
3、 根據權利要求1所述的基于重力場模型的管道狀物體拉直的圖像處理方 法，其特征是，所述的邊界條件是指在外壁三維有限元模型的一端取一橫截面， 設定位于該橫截面上的所有節點的位移和速度均為0。
4、 根據權利要求1所述的基于重力場模型的管道狀物體拉直的圖像處理方 法，其特征是，所述的重力場載荷是指設定垂直于施加邊界約束的橫截面，且指 向外壁三維有限元模型內側的方向為重力場方向。
5、 根據權利要求1所述的基于重力場模型的管道狀物體拉直的圖像處理方法，其特征是，所述的非線性有限元變形計算求解是指列出虛功方程，使用 Newton-R鄰hson方法、修正的Newton-Raphson方法或準Newton方法中的一種 計算虛功方程的極值點。
6、 根據權利要求1所述的基于重力場模型的管道狀物體拉直的圖像處理方 法，其特征是，所述的虛功方程是指變形體中任意滿足平衡的力系在任意滿足 協調條件的變形狀態上作的虛功等于零，即體系外力的虛功與內力的虛功之和等 于零。
7、 根據權利要求1所述的基于重力場模型的管道狀物體拉直的圖像處理方 法，其特征是，所述的非剛體配準方法是指薄盤樣條法、近似薄盤樣條法、體 積樣條法或彈性體樣條法中的一種。
8、 根據權利要求1所述的基于重力場模型的管道狀物體拉直的圖像處理方法，其特征是，所述的拉直變形后的外壁三維圖像數據是指利用獲取的拉直變形后的內壁三維圖像數據，通過膨脹腐蝕或非剛體配準等方法中的一種獲得的拉 直變形后的管道狀物體的內部空腔三維圖像數據、拉直變形后的管道狀物體的管 壁三維圖像數據的總和。
9、 根據權利要求1所述的基于重力場模型的管道狀物體拉直的圖像處理方法，其特征是，所述的長度失真校正處理，具體步驟包括1) 利用第七步得到的節點集合iV^A^A^L ,A^計算其中每一點 A^(A; = 1，2,L ,w-1)與下一相鄰節點iV^之間的直線段在拉直變形后的外壁三維 圖像數據的中心路徑上的投影長度《；2) 利用第七步得到的節點集合iV'^iV;,A^L ,《}計算其中每一點<formula>formula see original document page 4</formula>與下一相鄰節點《+1之間的直線段在外壁三維圖像數據的中 心路徑上的投影長度《；3)計算第l)步和第2)步所得的投影長度的差<formula>formula see original document page 4</formula>(7)作為中心路徑點(^+1所在數據層的長度失真的估計，使用基于特征點的非剛體配準方法進行長度失真校正對于長度失真校正后的內壁三維圖像數據的第 A:(A: = 1,2，L ")層，取點0c,乂A:)作為目標特征點，取點(x，;;,m)作為源特征點，其中<formula>formula see original document page 4</formula>附<formula>formula see original document page 4</formula>(8)if hi.其中x和_y為拉直變形后的外壁三維圖像數據的中心路徑在第A:層的坐標 值，根據源特征點與目標特征點的對應關系，求得每一對特征點之間的變換向量， 再對旋轉失真校正后的內壁三維圖像數據利用薄盤樣條插值，得到旋轉失真校正 后的內壁三維圖像數據中任意一點到長度失真校正后的內壁三維圖像數據中對 應點的映射關系，使用基于特征點的非剛體配準的方法實現長度失真校正，得到 長度失真校正后的內壁三維圖像數據。
10、根據權利要求1所述的基于重力場模型的管道狀物體拉直的圖像處理方 法，其特征是，所述的旋轉失真校正處理，具體步驟包括1) 對于管道狀物體三維數據中心路徑上的每一點C,，在外壁拉直變形后的有限元模型中選擇與之距離最小的節點iV,;2) 平移中心路徑，對于平移后中心路徑上的每一點C;，再次在外壁拉直變 形后的有限元模型中選擇與之距離最小的節點《；3) 在外壁拉直變形分析前模型中，可以分別找到與iV,和《一一對應的節點4) 用四組節點中序號相同的節點可以構建兩組向量，5 = $-g和uu1 uu 111 訓 uul4 uui分別計算^/,與下一相鄰向量[/,+1之間的角度《、《與下一相鄰向5) 以角度差z^,;《-《作為中心路徑點C,所在數據層的旋轉失真的估計，進行旋轉失真校正；6) 對拉直變形后的內壁三維圖像數據第二層以后的每一層，進行旋轉變換， 得到旋轉失真校正后的內壁三維圖像數據對拉直變形后的內壁三維圖像數據第 /(/ = 1,2,L ,")層的任意一點(x，力，與之對應的旋轉失真校正后的內壁三維圖像 數據點的坐標(rx，o;)滿足LL(x - cx) sin(z!D — (y _ cy) cos(zlD + cy」,if z' > 1,其中(cx,c力為拉直變形后的外壁三維圖像數據的中心路徑在第/層的坐標。if / = 1， if z' > 全文摘要
一種圖像處理技術領域的基于重力場模型的管道狀物體拉直的圖像處理方法，包括以下步驟對管道狀物體的內壁三維圖像數據和外壁三維圖像數據進行網格劃分并生成外壁三維有限元模型；對外壁三維有限元模型施加邊界條件和重力場載荷；進行非線性有限元拉直變形分析；使用基于特征點的非剛體配準的方法實現內壁三維圖像數據的拉直；采用中心路徑生成算法計算外壁三維圖像數據的中心路徑和拉直變形后的外壁三維圖像數據的中心路徑；進行旋轉失真校正處理和長度失真校正處理。本發明能夠進一步獲得優化的管道狀物體拉直后的圖像，避免了現有技術中基于電場模型和中心路徑拉直方法所造成的拉直算法失效問題。
文檔編號G06T15/00GK101615299SQ20091005560
公開日2009年12月30日 申請日期2009年7月30日 優先權日2009年7月30日
發明者雷 李, 俊 趙 申請人:上海交通大學