TWI594732B - Three-dimensional median filter applied to computed tomography - Google Patents

Three-dimensional median filter applied to computed tomography Download PDF

Info

Publication number
TWI594732B
TWI594732B TW104143969A TW104143969A TWI594732B TW I594732 B TWI594732 B TW I594732B TW 104143969 A TW104143969 A TW 104143969A TW 104143969 A TW104143969 A TW 104143969A TW I594732 B TWI594732 B TW I594732B
Authority
TW
Taiwan
Prior art keywords
noise
dimensional
median
image
coordinate
Prior art date
Application number
TW104143969A
Other languages
English (en)
Other versions
TW201722357A (zh
Inventor
Bing-Guo Weng
Ying-Yi Wu
Yi-Fu Tang
Lan-Rong Dong
Original Assignee
Nat Chung-Shan Inst Of Science And Tech
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nat Chung-Shan Inst Of Science And Tech filed Critical Nat Chung-Shan Inst Of Science And Tech
Priority to TW104143969A priority Critical patent/TWI594732B/zh
Publication of TW201722357A publication Critical patent/TW201722357A/zh
Application granted granted Critical
Publication of TWI594732B publication Critical patent/TWI594732B/zh

Links

Landscapes

  • Apparatus For Radiation Diagnosis (AREA)
  • Image Processing (AREA)

Description

應用於電腦斷層影像之三維中值濾波方法
本發明係關於一種電腦斷層影像的影像處理方法，特別是關於一種電腦斷層影像之三維中值濾波方法。
電腦斷層掃描Computer Tomography是用X光來產生影像，X光由X球管產生，當X光穿透不同密度的被檢測物時，X光會因組織密度的差異而產生不同程度的衰減，這些強弱程度不同的X光，會由閃爍體(CsI)吸收並轉換成可見光，可見光會被一個一維或二維的CCD或CMOS影像感測器吸收，最後得到被拍攝物體的X光投影影像，而X光管和偵測器則會持續在被檢測體的橫斷面上做旋轉，每轉一定角度就拍攝一張影像，直到完成旋轉180度以上的多方向影像拍攝；使用一維偵檢器時X球管需以螺旋方式對被檢測物旋轉多圈，使用二維偵檢器配合錐狀束X球管只需對被檢測物旋轉一圈即可，之後電腦會透過演算法將獲得的X光影像重建為二維的電腦斷層切面影像，二維切面影像堆疊起來可成為三維影像，三維影像係由二維切面影像合成而產生，因此往往會因為取像的精確度限制而摻雜斑點雜訊(impulsive noise)。
目前許多的影像處理技術提供多種方法解決雜訊問題，但如果目的只是要把雜訊去除，這樣的處理固然可以有效的在視覺上消除每張斷層切面影像(二維影像)的斑點雜訊，這些個別經影像處理的二維切面影像堆疊成為三維影像時，二維切面影像會造成三維影像成像時物體(objects)切片組合的不連續，這些不連續是肇因於影像處理後物體切片間的邊緣無法對齊，因此被檢測物3D重建後的三維影像會產生輪廓失真。
因此目前業界極需發展出一種三維電腦斷層影像的影像處理方法，可消除三維影像堆疊成像時物體(objects)切片組合的不連續，如此一來，方能同時兼具消除雜訊與避免三維電腦斷層影像建構時產生輪廓失真的問題。
鑒於上述悉知技術之缺點，本發明之主要目的在於提供一種應用於電腦斷層影像之三維中值濾波方法，整合一組可建構出三維影像的二維切面電腦斷層影像、一雜訊的分佈狀態、一濾波演算法及一中位值影像資料等，以建構出消除雜訊與避免輪廓失真的三維電腦斷層影像。
為了達到上述目的，根據本發明所提出之一方案，提供一種應用於電腦斷層影像之三維中值濾波方法，步驟包括：(A)讀取一組可建構出三維影像的二維切面電腦斷層影像；(B)標定該組二維切面電腦斷層影像所有雜訊的座標位 置；(C)分析每一點雜訊之其所在二維切面電腦斷層影像及相鄰的二維切面電腦斷層影像的雜訊分佈狀態；(D)根據雜訊的分佈狀態進行一濾波演算法而得一中位值影像資料，並以該中位值影像資料取代雜訊；(E)將該組二維切面電腦斷層影像的雜訊處理後，建構出三維影像。
步驟(B)中該雜訊的座標位置的標定，為給予位於第Z張二維切面電腦斷層影像的該雜訊一個笛卡兒座標值(x,y,z)，其中z係代表第Z張二維切面電腦斷層影像，x、y代表雜訊在第Z張二維切面電腦斷層影像平面的二維座標，藉由雜訊的笛卡兒座標值(x,y,z)，可精確定位出雜訊在該組可建構出三維影像的二維切面電腦斷層影像的位置。
步驟(C)中雜訊所在的二維切面電腦斷層影像及其相鄰的二維切面電腦斷層影像，指的是第Z張有雜訊的二維切面電腦斷層影像及第Z-1、Z+1張的二維切面電腦斷層影像；本發明分析處理的雜訊，不單純處理二維的雜訊資訊，而是分析雜訊在第Z-1、Z及Z+1張二維切面電腦斷層影像上的3維的分佈狀態。
根據雜訊的3維的分佈狀態，可將雜訊分為3種類型(但不以此為限)來做後續的影像處理，本發明提出三種區域範圍(但不以此為限)為第一等第笛卡兒座標區域、第二等第笛卡兒座標區域及第三等第笛卡兒座標區域；當雜訊點的3維的分佈狀態落入第一等第笛卡兒座標區域的範圍時，本發 明提出的濾波演算法為：I1(0,0,0)=median({I(x,y,x))|(x,y,z)S1})其中，S1表示第一等第笛卡兒座標區域，包含(0,0,0)、(1,0,0)、(0,1,0)、(0,0,1)、(-1,0,0)、(0,-1,0)、(0,0,-1)，其中，S1：第一等第笛卡兒座標區域、I1：執行濾波演算法後，雜訊在座標位置的中位值影像資料，(0,0,0)為已被標定的雜訊點在第一等第笛卡兒座標區域的座標；當雜訊的3維的分佈狀態落入第二等第笛卡兒座標區域的範圍時，本發明提出的濾波演算法為：I2(0,0,0)=median({I(x,y,x))|(x,y,z)S2})其中，S2表示第二等第笛卡兒座標區域，包含(0,0,0)、(1,0,0)、(0,1,0)、(0,0,1)、(-1,0,0)、(0,-1,0)、(0,0,-1)、(0,1,1)、(0,-1,1)、(0,1,-1)、(0,-1,-1)、(1,1,0)、(1,-1,0)、(-1,1,0)、(-1,-1,0)、(1,0,1)、(-1,0,1)、(1,0,-1)、(-1,0,-1)，其中，S2：第二等第笛卡兒座標區域、I2：執行濾波演算法後，雜訊在座標位置的中位值影像資料，(0,0,0)為已被標定的雜訊點在第二等第笛卡兒座標區域的座標；當雜訊的3維的分佈狀態落入第三等第笛卡兒座標區域的範圍時，本發明提出的濾波演算法為：I3(0,0,0)=median({I(x,y,x))|(x,y,z)S3})其中，S3表示第三等第笛卡兒座標區域，包含(0,0,0)、(1,0,0)、(0,1,0)、(0,0,1)、(-1,0,0)、(0,-1,0)、(0,0,-1)、(0,1,1)、(0,-1,1)、(0,1,-1)、(0,-1,-1)、(1,1,0)、(1,-1,0)、(-1,1,0)、(-1,-1,0)、(1,0,1)、(-1,0,1)、(1,0,-1)、(-1,0,-1)、(1,1,1)、(-1,1,1)、(1,-1,1)、(-1,-1,1)、(1,1,-1)、(1,-1,-1)、(-1,1,-1)、 (-1,-1,-1)，其中，S3：第三等第笛卡兒座標區域、I3：執行濾波演算法後，雜訊在座標位置的中位值影像資料，(0,0,0)為已被標定的雜訊點在第三等第笛卡兒座標區域的座標。
以上之概述與接下來的詳細說明及附圖，皆是為了能進一步說明本創作達到預定目的所採取的方式、手段及功效。而有關本創作的其他目的及優點，將在後續的說明及圖式中加以闡述。
S101-S105‧‧‧步驟
S1‧‧‧第一等第笛卡兒座標區域
S2‧‧‧第二等第笛卡兒座標區域
S3‧‧‧第三等第笛卡兒座標區域
第一圖係為本發明一種應用於電腦斷層影像之三維中值濾波方法流程圖；第二圖係為本發明第一等第笛卡兒座標區域示意圖；第三圖係為本發明第二等第笛卡兒座標區域示意圖；第四圖係為本發明第三等第笛卡兒座標區域示意圖。
以下係藉由特定的具體實例說明本創作之實施方式，熟悉此技藝之人士可由本說明書所揭示之內容輕易地了解本創作之優點及功效。
各種立體透視醫學影像(如：CT,MRI,PET等)往往會因為取像的精確度限制造成三維的合成影像會摻雜斑點雜 訊(impulsive noise)，本發明消除斑點雜訊的方法是採用中值濾波法對每一張二維電腦斷層影像進行雜訊濾除，但若單純使用中值濾波法對每一張二維電腦斷層影像進行雜訊濾除，固然可以在視覺上消除每張斷層影像的斑點雜訊，但是卻會造成三維影像成像時物體(objects)切片組合的不連續，使得濾波後物體切片間的邊緣無法對齊。
請參閱第一圖，為本發明一種應用於電腦斷層影像之三維中值濾波方法流程圖。如圖所示，本發明所提供一種應用於電腦斷層影像之三維中值濾波方法，步驟包括：(A)讀取一組可建構出三維影像的二維切面電腦斷層影像S101；(B)標定該組二維切面電腦斷層影像所有雜訊的座標位置S102；分析每一點雜訊之其所在二維切面電腦斷層影像及相鄰的二維切面電腦斷層影像的分佈狀態S103；根據雜訊的分佈狀態進行一濾波演算法而得一中位值影像資料，並以該中位值影像資料取代雜訊S104；將該組二維切面電腦斷層影像的雜訊處理後，建構出三維影像S105。
本發明的濾波演算法是利用中值濾波器進行雜訊濾除，中值濾波器(median filter)在脈衝雜訊(impulse noises)出現時特別有用，因為脈衝雜訊看起來，像是疊加在影像上的白點和黑點，所以又稱為胡椒鹽雜訊(salt-and-pepper noise)；本發明利用中值濾波器的濾波演算法，係為使用遮罩之排序運算，遮罩之排序運算是一種非 線性的濾波方式，首先對在遮罩裡所有像素的灰階值大小做排序，接著選取排序在中間的灰階值(或稱中值灰階值、中位值灰階值)的影像資料，再將此中位值(灰階)影像資料代入所對應的像素(雜訊影像資料)中，如此便完成一個像素的濾波動作，此法不但可以濾掉影像中突起的高頻雜訊部份，對於影像的邊緣也能夠給予適當的保留，再將像素的值(雜訊的影像資料)用該像素近鄰灰階的中位值影像資料(經濾波演算法運算)來取代，本發明即是根據雜訊大小，選擇不同遮罩的中值濾波器，遮罩的形式是三個等第的笛卡兒座標區域。
請參閱第二圖，為本發明第一等第笛卡兒座標區域示意圖、請參閱第三圖，為本發明第二等第笛卡兒座標區域示意圖、請參閱第四圖，為本發明第三等第笛卡兒座標區域示意圖。本發明實施例步驟如下：
步驟(1)首先讀取一組可建構出三維影像的二維切面電腦斷層影像，以(x,y,z)代表第z張二維切面影像的(x,y)座標值，將二維切面影像的平面數位化，以方便位置標定及影像資料的數位化處理。
步驟(2)定義出三種不同區域大小的笛卡兒座標區域，，包含第一等第笛卡兒座標區域、第二等第笛卡兒座標區域、第三等第笛卡兒座標區域，這三種不同區域大小的笛卡兒座標區域裡的座標不同於步驟(1)的座標(x,y,z)，其中，第 一等第笛卡兒座標區域(S1)包括(0,0,0)、(1,0,0)、(0,1,0)、(0,0,1)、(-1,0,0)、(0,-1,0)、(0,0,-1)(如圖二所示)，第二等第笛卡兒座標區域(S2)包括(0,0,0)、(1,0,0)、(0,1,0)、(0,0,1)、(-1,0,0)、(0,-1,0)、(0,0,-1)、(0,1,1)、(0,-1,1)、(0,1,-1)、(0,-1,-1)、(1,1,0)、(1,-1,0)、(-1,1,0)、(-1,-1,0)、(1,0,1)、(-1,0,1)、(1,0,-1)、(-1,0,-1)(如圖三所示)，第三等第笛卡兒座標區域(S3)包括(0,0,0)、(1,0,0)、(0,1,0)、(0,0,1)、(-1,0,0)、(0,-1,0)、(0,0,-1)、(0,1,1)、(0,-1,1)、(0,1,-1)、(0,-1,-1)、(1,1,0)、(1,-1,0)、(-1,1,0)、(-1,-1,0)、(1,0,1)、(-1,0,1)、(1,0,-1)、(-1,0,-1)、(1,1,1)、(-1,1,1)、(1,-1,1)、(-1,-1,1)、(1,1,-1)、(1,-1,-1)、(-1,1,-1)、(-1,-1,-1)(如圖四所示)。
步驟(3)檢查第z片與前後片(第z+1張、第z-1張)共三張二維切面斷層影像之雜訊分佈，並記錄下所有雜訊座標位置的座標點(x,y,z)，其中，雜訊通常不會是一個個體，一個雜訊的分佈可能包含許多雜訊顆粒，因此雜訊座標位置的座標點(x,y,z)可利用質心點(不以此為限)的計算方式，計算出某一群落雜訊顆粒的質心點，以該質心點座標為雜訊座標位置的座標點(x,y,z)，以方便標定每一點雜訊在該組可建構出三維影像的二維切面電腦斷層影像的位置，但在進行後續濾波演算法(中值濾波)時，每一雜訊點在其所屬的笛卡兒座標區域中的座標皆設置為(0,0,0)。
步驟(4)依據雜訊分佈選擇中值濾波器的遮罩大小與相鄰點位置，當雜訊分佈散落於第一等第笛卡兒座標區域時，中值濾波器的遮罩大小為7，而相鄰點位置(雜訊分佈區域)的選擇為S1；當雜訊散落於第二等第笛卡兒座標區域時，中值濾波器的遮罩大小為19，而相鄰點位置(雜訊分佈區域)的選擇為S2；當雜訊散落於第三等第笛卡兒座標區域時，中值濾波器的遮罩大小為27，而相鄰點位置(雜訊分佈區域)的選擇為S3。
步驟(5)依雜訊分佈的區域範圍座落於哪一個等第笛卡兒座標區域，執行不同的濾波演算法(中值濾波)，若雜訊分佈位於第一等第笛卡兒座標區域時，使用式(一)完成濾波：I1(0,0,0)=median({I(x,y,x))|(x,y,z)S1})...式(一)式(一)中，I1為執行濾波演算法後，雜訊在座標位置的中間灰階值(中位值影像資料)，座標(0,0,0)為雜訊座標位置在第一等第笛卡兒座標區域(S1)的座標，其餘座標為其他雜訊顆粒(雜訊分布)在第一等第笛卡兒座標區域(S1)的座標；若雜訊分佈位於第一等第笛卡兒座標區域時，使用式(二)完成濾波：I2(0,0,0)=median({I(x,y,x))|(x,y,z)S2})...式(二)式(二)中，I2為執行濾波演算法後，雜訊在座標位置的中間灰階值(中位值影像資料)，座標(0,0,0)為雜訊座標位置在第二等第笛卡兒座標區域(S2)的座標，其餘座標為其他雜訊顆 粒(雜訊分布)在第二等第笛卡兒座標區域(S2)的座標；若雜訊分佈位於第三等第笛卡兒座標區域時，使用式(三)完成濾波：I3(0,0,0)=median({I(x,y,x))|(x,y,z)S3})...式(三)式(三)中，I3為執行濾波演算法後，雜訊在座標位置的中間灰階值(中位值影像資料)，座標(0,0,0)為雜訊座標位置在第二等第笛卡兒座標區域(S3)的座標，其餘座標為其他雜訊顆粒(雜訊分布)在第三等第笛卡兒座標區域(S3)的座標。
步驟(6)重複步驟(4)-(5)，直到所有雜訊處理完畢。
步驟(7)將所有的二維切面電腦斷層影像的雜訊處理後，依順序建構出三維影像。
上述之實施例僅為例示性說明本創作之特點及功效，非用以限制本創作之實質技術內容的範圍。任何熟悉此技藝之人士均可在不違背創作之精神及範疇下，對上述實施例進行修飾與變化。因此，本創作之權利保護範圍，應如後述之申請專利範圍所列。
S101-S105‧‧‧步驟

Claims (7)

  1. 一種應用於電腦斷層影像之三維中值濾波方法，步驟包括：(A)讀取一組可建構出三維影像的二維切面電腦斷層影像；(B)標定該組二維切面電腦斷層影像所有雜訊的座標位置；(C)分析每一點雜訊之其所在二維切面電腦斷層影像及相鄰的二維切面電腦斷層影像的分佈狀態；(D)根據雜訊的分佈狀態進行一濾波演算法而得一中位值影像資料，並以該中位值影像資料取代雜訊；(E)不分析雜訊以外其他位置影像資料，將該組二維切面電腦斷層影像的雜訊處理後，建構出三維影像；其中，該雜訊係為複數雜訊顆粒。
  2. 如申請專利範圍第1項所述之應用於電腦斷層影像之三維中值濾波方法，其中，步驟(B)中該雜訊的座標位置之座標係依位於第Z張二維切面電腦斷層影像的該雜訊給予一笛卡兒座標值(x,y,z)。
  3. 如申請專利範圍第2項所述之應用於電腦斷層影像之三維中值濾波方法，其中，該雜訊的分佈狀態係為該雜訊在第Z-1、Z及Z+1張二維切面電腦斷層影像上的分佈狀態。
  4. 如申請專利範圍第2項所述之應用於電腦斷層影像之三維中值濾波方法，其中，步驟(C)中該雜訊的分佈狀態係依一 第一等第笛卡兒座標區域、一第二等第笛卡兒座標區域、一第三等第笛卡兒座標區域來分類。
  5. 如申請專利範圍第4項所述之應用於電腦斷層影像之三維中值濾波方法，其中，步驟(D)中屬於該第一等第笛卡兒座標區域的雜訊分佈進行的該濾波演算法係為：I1=median({I(x,y,x))|(x,y,z)S1})S1：第一等第笛卡兒座標區域、I1：執行濾波演算法後，雜訊在座標位置的中位值影像資料。
  6. 如申請專利範圍第4項所述之應用於電腦斷層影像之三維中值濾波方法，其中，步驟(D)中屬於該第二等第笛卡兒座標的雜訊分佈進行的該濾波演算法係為：I2=median({I(x,y,x))|(x,y,z)S2})S2：第二等第笛卡兒座標區域、I2：執行濾波演算法後，雜訊在座標位置的中位值影像資料。
  7. 如申請專利範圍第4項所述之應用於電腦斷層影像之三維中值濾波方法，其中，步驟(D)中屬於該第三等第笛卡兒座標的雜訊分佈進行的該濾波演算法係為：I3=median({I(x,y,x))|(x,y,z)S3})S3：第三等第笛卡兒座標區域、I3：執行濾波演算法後，雜訊在座標位置的中位值影像資料。
TW104143969A 2015-12-28 2015-12-28 Three-dimensional median filter applied to computed tomography TWI594732B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
TW104143969A TWI594732B (zh) 2015-12-28 2015-12-28 Three-dimensional median filter applied to computed tomography

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
TW104143969A TWI594732B (zh) 2015-12-28 2015-12-28 Three-dimensional median filter applied to computed tomography

Publications (2)

Publication Number Publication Date
TW201722357A TW201722357A (zh) 2017-07-01
TWI594732B true TWI594732B (zh) 2017-08-11

Family

ID=60047991

Family Applications (1)

Application Number Title Priority Date Filing Date
TW104143969A TWI594732B (zh) 2015-12-28 2015-12-28 Three-dimensional median filter applied to computed tomography

Country Status (1)

Country Link
TW (1) TWI594732B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108295384B (zh) 2017-01-11 2020-02-28 南京中硼联康医疗科技有限公司 基于医学影像的组织元素质量比例解构方法及几何模型建立方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040081340A1 (en) * 2002-10-28 2004-04-29 Kabushiki Kaisha Toshiba Image processing apparatus and ultrasound diagnosis apparatus
TW201109949A (en) * 2009-09-01 2011-03-16 Univ Nat Pingtung Sci & Tech Density-based data clustering method

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040081340A1 (en) * 2002-10-28 2004-04-29 Kabushiki Kaisha Toshiba Image processing apparatus and ultrasound diagnosis apparatus
TW201109949A (en) * 2009-09-01 2011-03-16 Univ Nat Pingtung Sci & Tech Density-based data clustering method

Also Published As

Publication number Publication date
TW201722357A (zh) 2017-07-01

Similar Documents

Publication Publication Date Title
US10274439B2 (en) System and method for spectral x-ray imaging
JP5858625B2 (ja) 補正版再構成コンピュータ断層像生成方法
JP4354484B2 (ja) 内部構造の撮像
CN103679642B (zh) 一种ct图像金属伪影校正方法、装置及ct设备
KR101531440B1 (ko) 표면 스캔 정보를 이용한 3차원 치과용 x 레이 데이터 세트로부터의 아티팩트의 감소 및 제거
JP4175468B2 (ja) 画像処理方法、画像処理プログラム及びコンピュータ読取可能な記録媒体
JP2007050259A (ja) ボリュームデータ再構成後の断層撮影３ｄ画像のフィルタリング方法
ES2706749T3 (es) Método para procesar datos de imagen que representan un volumen tridimensional
JP6273291B2 (ja) 画像処理装置および方法
JP2007530132A (ja) コンピュータ断層撮影画像を向上させるための方法、コンピュータプログラムプロダクトおよび装置
JP2008043758A (ja) オパシティの検出のために放射線画像を処理するための方法
KR20130128124A (ko) 파노라마 영상 데이터 제공 방법 및 장치
JP2016145778A (ja) Ｘ線検査装置及びｘ線検査方法
US9672641B2 (en) Method, apparatus, and computer readable medium for removing unwanted objects from a tomogram
CN107106107B (zh) 用于处理断层合成数据的方法和系统
JP7078642B2 (ja) 物体の複数の構成要素の３ｄモデルデータを取得する方法
Körner et al. Increasing throughput in x-ray computed tomography measurement of surface topography using sinogram interpolation
Mouton et al. A novel intensity limiting approach to metal artefact reduction in 3D CT baggage imagery
CN106108932B (zh) 全自动肾脏感兴趣区提取装置与方法
Johari et al. Metal artifact suppression in dental cone beam computed tomography images using image processing techniques
US9867586B2 (en) Stereo X-ray tube based suppression of outside body high contrast objects
TWI594732B (zh) Three-dimensional median filter applied to computed tomography
JP6525837B2 (ja) 製品の欠陥検出方法
CN109870471B (zh) 一种单光栅侦测的锥束ct角度序列散射获取方法
KR20170109519A (ko) 핵의학 화상 분석 기술