CT斷層圖像重建算法研究

年4月19日CT斷層圖像重建算法研究文檔僅供參考CT斷層圖像重建算法研究專業(yè)：通信工程姓名：劉明帥指導教師：駱巖紅摘要CT技術是一種融合了射線光電子學、信息學、微電子學等學科的新興技術，因為其先進的無損檢測技術，因此被廣泛地應用于醫(yī)學、航天、生物等多個領域。隨著科技的進步，圖像重建技術開始應用于X射線中，這是數(shù)字圖像處理的一個重大進步。如何能重建出高質量的圖像，取決于所采用的重建算法。從圖像重建的角度來看，主要分為解析法與迭代法。解析法是利用解析、變換重建公式來構建重建圖像。它具有容易實現(xiàn),速度較快,且能重建出高質量的圖像的特點,可是對投影數(shù)據(jù)完備性要求高。迭代法是利用求解線性方程組來重建圖像，它能夠在投影數(shù)據(jù)信噪較低條件下，獲得高質量圖像。本文將從原理、應用、與優(yōu)缺點的角度來分析兩種算法，重點對解析法中的濾波反投影算法從平行束與扇束投影方式進行研究，最后經過VisualC++與MATLAB軟件相結合的方式對圖像重建，并分析各參數(shù)對重建圖像的影響。關鍵字：CT技術圖像重建算法濾波反投影算法AbstractCTtechnologyisaemergingtechnologythatblendoftheRayoptoelectronics,microelectronicsandinformaticssubject.Becauseofitsadvancednondestructivetestingtechnology,itiswidelyusedinmedical,aerospace,biologicalandotherfields.Withtheprogressofscienceandtechnology,ImagereconstructiontechnologyisappliedtotheXray,Thisisamajorprogressofdigitalimageprocessing.Howtorebuildthehighqualityimages,dependsonthereconstructionalgorithmyouadopt.Fromtheperspectiveofimagereconstruction,

itmainlydividedintotheanalyticalmethodanditerationmethod.Analyticalmethodusetheanalysisandtransformformulatobuildimagereconstruction.Ithasthecharacteristicsofimplementatingeasilyandfast,andreconstructingouthighqualityimages,butthedemandoftheprojectiondataishigh.Iterativemethodisusedtosolvethelinearsystemofequationstoreconstructionimage,theprojectiondataundertheconditionoflowsignal-to-noise,itcangethighqualityimage.Thisarticlewewillbefromthepointofviewoftheprinciple,application,andtheadvantagesanddisadvantagestoanalysisthetwokindsofalgorithms,

focusingonstudyingtheanalyticalmethodoffilterbackprojectionalgorithmfromtheparallelbeamandfanbeamprojectionmethods,finally,combiningthesoftwareofVisualc++withMATLABsoftwaretoimagereconstruction,andanalyzestheinfluenceofvariousparametersonthereconstructionimageKeywords:CTtechnologyimagereconstructionalgorithmFilteredBackprojectionAlgorithm目錄第一章緒論 -5-1.1CT技術與圖像重建概述 -6-1.2CT和重建技術的發(fā)展及研究現(xiàn)狀 -6-1.3研究的目的與意義 -9-第二章CT成像原理和圖像重建算法 -9-2.1CT成像原理與系統(tǒng)組成 -9-2.2CT成像系統(tǒng)掃描方式的發(fā)展 -10-2.3CT斷層圖像原理 -12-2.4圖像重建算法概述 -13-2.4.1解析類方法 -13-2.4.2傳統(tǒng)迭代類方法 -14-第三章CT圖像重建算法實現(xiàn)原理的研究 -14-3.1圖像重建系統(tǒng)中的數(shù)學概念及變換 -14-3.1.1投影與反投影 -14-3.1.2Radon變換及其反變換 -15-3.1.3傅里葉變換 -16-3.1.4中心切片定理 -17-3.2解析類重建算法 -18-3.2.1直接傅里葉變換算法 -18-3.2.2反投影重建算法 -19-3.3迭代類重建算法 -21-3.3.1代數(shù)迭代重建算法 -22-(1)ART算法 -22-(2)同時代數(shù)重建算法 -22-3.3.2影響代數(shù)迭代重建算法的因素 -23-3.3.3ART重建算法與SART -26-3.3.4統(tǒng)計迭代重建算法 -27-（1）最小二乘圖像重建算法 -27-第四章濾波反投影重建算法 -30-4.1平行束濾波反投影重建算法 -31-4.1.1卷積反投影重建算法 -31-4.1.2濾波函數(shù)的選擇 -32-4.2扇形束濾波反投影算法 -34-4.2.1等夾角扇形束濾波反投影算法 -35-4.2.2等間距扇形濾波反投影算法 -36-4.3濾波反投影重建算法的軟件實現(xiàn) -37-4.4濾波反投影算法和ART算法的對比 -49-第五章總結與展望 -50-5.1總結 -50-5.2新興的迭代算法 -51-5.3展望 -52-參考文獻 -53-第一章緒論1.1CT技術與圖像重建概述所謂的CT技術就是X射線計算機斷層成像，它是一種新興技術，發(fā)展于20世紀80年代，它將射線光電子學、信息科學、微電子學、計算機科學等學科結合在一起。我們知道，當X射線照射不同物體時，每個物體對這種射線的吸收和透射率不同，而重建正是利用這種原理，射線照射后，利用探測器進行接收，這樣我們就能夠根據(jù)衰減數(shù)得到其分布圖像，這就是CT成像技術的基礎。CT技術在對物體進行檢測時，不用破壞物體內部結構,正是由于這種無損檢測技術，因此被廣泛應用于醫(yī)學、生物、航天、航空等多個領域。圖像重建是由物體的截斷面向該平面做投影，根據(jù)投影所得的函數(shù)來重建截斷面的過程。隨著時間的推移，人類在科學技術上也有了重大進步，特別是計算機技術的高速發(fā)展，也推動了圖像重建的發(fā)展，在醫(yī)學領域應用最為顯著，它大大的豐富了對于人體內臟器官檢查的手段，為更加準確地診斷疾病提供了強有力的依據(jù)。根據(jù)原始數(shù)據(jù)獲取方法和重建原理不同，可分為透射斷層重建成像、發(fā)射斷層重建成像、反射斷層重建成像[1]。1.2CT和重建技術的發(fā)展及研究現(xiàn)狀早在1895年，倫琴發(fā)現(xiàn)了X射線，這就是CT技術發(fā)展的早期萌芽。自此之后，人們也意識到CT技術在成像上有很好地發(fā)展前景。果如其然，不就之后人們就看到了諸如旋轉陽極X射線管、影像增強管等東西。然而，圖像重建技術應用于X放射線醫(yī)學卻是一個重大突破，對身體內臟疾病診斷有著非凡的意義。19，奧地利數(shù)學家J.Radon發(fā)表了圖像投影理論，可是由于當時科學技術受限，因此沒能將其實現(xiàn)。到了二十世紀六十年代，由于計算機技術的迅猛發(fā)展，圖像重建技術被不少學者用來創(chuàng)造性地探索研究，在這其間，英國物理學家Cormack提出了一種方案使得X射線投影在醫(yī)學領域中的應用變成可能。1967年，在英國EMI實驗中心從事計算機和重建技術研究的Hounsfield，她發(fā)現(xiàn)X射線從不同方向透過物體后，對這些衰減量做測試就能得到物體內部結構，經過不懈的研究，終于在1971年，她研制出了世界上第一臺臨床使用的計算機斷層掃描裝置，即CT機。1972年，這臺CT成功為一名婦女診斷。顱腦CT的安裝與使用，標志CT時代的到來。圖1-1和1-2分別表示CT掃描成像設備與臨床圖像。圖1-1CT掃描成像設備圖l-2CT臨床圖像20世紀80年代，能夠不間斷采集投影數(shù)據(jù)的螺旋式CT的誕生與使用，能獲得高質量的三維圖像。這期間，中國分別于1983年和1985年研制出了第一臺顱腦掃描裝置和第二代CT設備。1989年，中國又成功研制開發(fā)出了第一臺全身常規(guī)CT-。1992年，中國生產出了第一臺螺旋掃描設備。隨著人類的研究探索，多層螺旋在單層的基礎上發(fā)展起來，它具有覆蓋面積更廣，掃描速度更快，得到更高質量三維圖像，輻射低等特點。然而，被業(yè)界譽為CT技術的一重大突破的雙能CT于11月誕生了，它能以較低的劑量重建出具有很高時間分辨率的圖像。1.3研究的目的與意義圖像重建是利用物體截斷面的投影來重建它本身的圖像，它不需要將物體解剖，而正是由于這種先進的無損檢測技術，因此被廣泛地應用于檢測和觀察中?，F(xiàn)實生活中，有些物體在構建其本身圖像時不能破壞它的物理性，因此圖像重建特別在醫(yī)學領域中，它作為一種先進的檢測技術，能更準確地檢測人體內部的器官。圖像重建其實能夠看做是一類特殊的圖像復原技術?，F(xiàn)實中，有多種重建算法，不同的重建算法得出不同質量的重建圖像。當前，濾波反投影算法在CT斷層圖像重建中應用最為廣泛，因為它具有很高的精度，能夠快速實現(xiàn)，且它的基本算法容易在軟件和硬件上實現(xiàn)。第二章CT成像原理和圖像重建算法2.1CT成像原理與系統(tǒng)組成一般來說，CT成像過程大致可歸納為:用儀器掃描需要檢測的物體,由探測器接收透過物體的X射線,經過模數(shù)轉換后成為數(shù)字信號,即原始數(shù)據(jù)。然后經過射束硬化、解除零點漂移等預處理,以便獲得更精確投影數(shù)據(jù)；再經過圖像重建與處理后,可得到斷層圖像輸出顯示或存儲至設備。簡單地說，CT成像系統(tǒng)組成如下：(1)掃描設備，包括滑環(huán)和(平板)檢測器，X射線管，模數(shù)信號轉換系統(tǒng)等[2]；(2)原始數(shù)據(jù)處理與重建設備，即電子計算機系統(tǒng)；(3)圖像的存儲設備以及用來顯示的設備。CT成像過程以及其系統(tǒng)組成如圖2-1所示。2-1CT系統(tǒng)組成與成像過程2.2CT成像系統(tǒng)掃描方式的發(fā)展自1971年9月世界上第一臺CT裝置誕生以后，CT裝置很快推廣到其它領域，由于各領域的不同需要，系統(tǒng)的掃描方式也發(fā)生很大變化。當前為止，通用的掃描方式有四種：第一代CT掃描方式最簡單，就是平行束掃描。它由單一的射線管和在測試區(qū)相對應的探測器。掃描中，射線源和探測器同步旋轉和平移運動時，探測器記錄X射線衰減情況，在不同角度下，會產生不同投影數(shù)據(jù)，但該光片不能得出被檢測物體的衰減分布情況，需要分度旋轉和重復平移。如圖2-2所示。2-2第一代CT系統(tǒng)掃描方式2-3第二代CT系統(tǒng)掃描方式第二代與第一代CT掃描方式不同的是平移的次數(shù)減少了。該系統(tǒng)包括一個產生窄角扇束射線源和小型的探測器陣列，由于扇束的扇角比較小，因此需要將射線源和探測器旋轉前進行一定筆束平移。如圖2-3所示。第三代CT系統(tǒng)使用扇束掃描方式，它是使用產生大角度的X射線源，且使用探測器數(shù)量較多的陣列。這樣，射線源與探測器組成的扇形能夠覆蓋整個被檢測物體，不需平移，只需旋轉就能夠。如圖2-4所示。第四代CT系統(tǒng)同樣也是扇形束掃描方式，探測器陣列使用固定環(huán)形探測器陣列，如圖2-5所示。它僅僅是射線源進行旋轉。被檢測物體放在環(huán)形探測器的中心，扇形束扇角的大小能夠根據(jù)檢測區(qū)域來適當調整，射線源每次進行一定分度的旋轉，最后獲得能夠重建被檢測物體內部圖像的投影數(shù)據(jù)。2-4第三代CT系統(tǒng)掃描方式2-5第四代CT系統(tǒng)掃描方式2.3CT斷層圖像原理根據(jù)投影圖像的重建原理以及方法，能夠分為不同的類型。我們主要研究透射成像原理，即斷層成像技術。斷層攝影成像的方式一般為發(fā)射接收，即用平行的X光線從不同的方向對物體進行照射，然后記錄每個方向的透射場，如圖2-6所示。2-6斷層圖像獲取示意圖2-7圖像在角下投影示意圖其中，X射線源產生平行的X射線對物體照射，設入射光強度為，接收器陣列得到的透射光度為。依次類推，X源和X射線接收器沿中心轉一個角度到，這樣就能夠一組投影數(shù)據(jù)為，其中的范圍是到，每隔一定間隔給予改變，這樣我們能夠得到投影向量。在數(shù)學上，圖像的投影能夠描述如下。如圖2-7所示，圖像函數(shù)為，穿過該線的一條線稱為射線。如果將函數(shù)在某條射線上的積分值化成一個集合，這個集合所表示的含義就是投影值。從原點向射線作一條垂線，此垂線作為新坐標的一個軸，并構成一個新坐標系，這樣能夠得出坐標系僅是坐標系旋轉角的結果，二者存在下列變換關系[3]： (2-1)因而射線積分表示為： == (2-2)斷層成像技術就是從不同射線角度，不同檢測器的位置的許多投影值重建原始圖像的過程，反應了處的密度。2.4圖像重建算法概述圖像重建在CT技術中發(fā)揮著重要作用。本質上說，它是按照采集后的數(shù)據(jù)，利用電子計算機來求解圖像矩陣中的像素，然后重新構建圖像的過程。從CT圖像重建的角度看,主要分為兩大類方法，即解析類方法與迭代類方法。2.4.1解析類方法解析類重建方法是利用解析、變換重建公式來構建重建圖像。當前，常見的解析類重建方法主要有:(1)濾波反投影算法(FBP)；(2)直接傅立葉變換算法(FBP)；(3)卷積反投影算法(CBP)。解析類算法因為其容易實現(xiàn),速度較快,且能重建出高質量的圖像的特點,因此CT系統(tǒng)被較為廣泛地應用起來，特別是在醫(yī)用CT系統(tǒng)中。然而濾波反投影算法應用最為廣泛?？墒?，其最大的不足之處在于對投影數(shù)據(jù)完備性的要求偏高,如果能夠在投影數(shù)據(jù)輸入給解析法之前，把不利于投影數(shù)據(jù)精確性的因素給與糾正，這樣就能夠得到滿意的重建圖像。2.4.2傳統(tǒng)迭代類方法迭代方法于1970年第一次被引入圖像重建的領域,它是利用數(shù)學級數(shù)迭代的原理來完成圖像重建的技術。迭代重建算法分別為代數(shù)迭代算法與統(tǒng)計迭代算法。代數(shù)迭代算法可分為:ART和SIRT型。統(tǒng)計迭代重建算法(SIR)是對投影數(shù)據(jù)進行準確性分析,經過精細的統(tǒng)計估計,逐步提高重建圖像的質量。迭代圖像重建技術的最大優(yōu)勢,在于它能夠在投影數(shù)據(jù)信噪比很低的情況下，甚至不完全的條件下，依然能夠獲得較高質量的重建圖像。然而迭代算法是逼近原始圖像，因此在實現(xiàn)時，運行時間相對較長且數(shù)據(jù)存儲量大。第三章CT圖像重建算法實現(xiàn)原理的研究3.1圖像重建系統(tǒng)中的數(shù)學概念及變換3.1.1投影與反投影我們知道，射線穿過物體之后，由于物體的吸收或散射作用，致使我們在檢測時會發(fā)現(xiàn)射線的強度會發(fā)生衰減。一般我們用衰減系數(shù)表示衰減程度。設一物質是分均勻的，一個面上衰減系數(shù)為，射線穿過該物質后，入射強度由變?yōu)?射線在面內的路徑如下圖3-1。3-1射線的衰減由Beer定律得出關系式： (3-1)在CT系統(tǒng)中，我們能夠把看作是射線在空氣中掃描時探測器測得的數(shù)據(jù)，是射線經過物體后衰減后數(shù)據(jù)，由于它們都是測量值，能夠得：(3-2)這樣我們就把上面的積分集合稱之為投影數(shù)據(jù)。將其推廣得出，如果射線從不同方向照射時，它所對應的路徑上投影數(shù)據(jù)都能得到，這樣構成一個投影數(shù)據(jù)集合。圖像重建就是利用投影數(shù)據(jù)集合來計算的過程。反投影的定義取決于投影的定義。但它不是投影運算的逆運算。數(shù)學語言來說，反投影算子不是投影算子的逆算子。3.1.2Radon變換及其反變換19，Radon變換的提出使CT的逆問題得到了解決，它也是CT圖像重建的理論基礎。如圖3-2所示。3-2Radon變換示意圖設直線L的方程為： (3-3)其中為直線到源點的距離，為軸正方向與直線垂線的角度，則與平面上直線L互相對應。記為要重建的圖像函數(shù)。有變換式如下， (3-4)稱該式為的Radon變換。其實這相當于一個算子，記為，它將空間中函數(shù)與空間中函數(shù)聯(lián)系在一起。因為是沿的積分，能夠認為空間上的一個點相應于空間上的一條直線L，用向量記，或者 (3-5)經過變換，上式變?yōu)? (3-6)這樣，相應的重建過程就是Radon反變換。3.1.3傅里葉變換一元連續(xù)函數(shù)傅里葉變換式定義為： (3-7)其中，。若給定的為連續(xù)函數(shù)，則反函數(shù)也存在，用傅里葉逆變換得到其原函數(shù)，即： (3-8)上倆式構成傅里葉變換對。上述式子推廣到二維中，設二元連續(xù)函數(shù)為，其傅里葉變換及其反變換定義式下： (3-9) (3-10)對于傅里葉變換，有一個重要定理，如下： (3-11)以上式子就是我們常見的卷積定理。3.1.4中心切片定理中心切片定理是斷層成像的理論基礎。二維圖像的中心切片定理指出：二維函數(shù)圖像在視角為時的投影的一維傅里葉變換,的二維傅里葉變換與探測器平行的方向，而且過原點的的一個切片。是切片與軸所成的夾角。原理如下圖3-3所示。3-3二維中心切片原理圖若探測器繞物體旋轉至少180度，物體的二維傅里葉變換所對應的探測器方向的中心片段就能夠將整個傅里葉空間進行覆蓋[4]。3.2解析類重建算法 解析類算法因其利用變換、重建公式，速度較快，實現(xiàn)容易，因此被廣泛應用。3.2.1直接傅里葉變換算法直接傅里葉變換算法是利用二維傅里葉變換來連接極變量函數(shù)與，它是直接利用投影定理。如下圖3-4。3-4投影參數(shù)說明圖將極坐標形式Radon變換改寫成直角坐標系下的形式，即： (3-12)其中表示平行束投影。進一步轉化有： (3-13)使用該方法進行圖像重建的過程是，先采集投影，再傅里葉變換，將投影放在傅里葉空間上，處理完所有投影后，在進行二維傅里葉逆變換，最后得到最終結果。但該算法一般不用于重建，僅作為理論研究，因為它的逆變換要得到星狀分布的變換插值為均勻網格形式，這是很困難的，同時它的計算量大，不能即刻成像，且重建后圖像質量比較差。3.2.2反投影重建算法反投影算法是最基本、最簡單的算法，它的基本思想是斷層內該點的密度值就是經過該平面上的一點的射線投影和。我們將“取投影”“反投影重建”“重建后圖像”看作一個系統(tǒng)，得出系統(tǒng)原理圖如下圖3-5所示。取投影取投影反投影重建原理重建后圖像3-5反投影系統(tǒng)原理圖設一位于坐標原點上的點源為斷層圖像上的唯一像點。當系統(tǒng)掃描方式為平移旋轉時，即射線先平移，再旋轉一定角度，知道累計轉角為為止。其中，為原坐標系與旋轉坐標系的夾角，這樣，投影位置可由來確定。如下圖3-6所示。3-6平移/旋轉掃描方式所用系統(tǒng)為固定坐標，為旋轉坐標系，為極坐標設為離散取值，當時，相應投影為：(3-14)若，相應投影為： (3-15)根據(jù)反投影公式的定義，點的圖像在坐標系表示為： (3-16)式中，，是投影數(shù)。上式反映的物理意義是：把經過某點所有的投影值做平均之后就是該點的密度值。如果有限區(qū)間將射線擴增至無限條，這樣就得到連續(xù)投影。這樣我們能夠得到投影表示式為： (3-17)式中，積分區(qū)間為，因為忽略射線硬化條件下，它與內投影值等效。在輸入圖像為點源條件下，我們得 (3-18)由式看出，反投影重建算法的系統(tǒng)擴展函數(shù)不是函數(shù)，因此系統(tǒng)不完美，圖像密度為零的點，重建后不一定為零，可能致圖像失真[5]。對于最為常見的濾波反投影算法，我將著重在下一章介紹。3.3迭代類重建算法由于計算機運算能力越來越強，使得迭代算法在應用中也越來越多地得到重視。現(xiàn)實中，迭代法重建圖像效果好，空間分辨率高，特別是當采集的數(shù)據(jù)不完全時，這一點體現(xiàn)更為突出。然而因為它是漸漸靠近原始圖像的，因此它在運算時間上會比較長，同樣帶來的是數(shù)據(jù)存儲量大。迭代法不是找到一個準確的解析式，這是與解析法最大的不同，一般來說，迭代算法的步驟是：先給斷層圖像賦一個初始估計值，根據(jù)此值算出理論投影值，將理論投影值與實際的相比較，按照一定的原則對原始圖像進行修正之后與理論值比較，在修正，如此循環(huán)，直到達到滿意的效果。其實概括地講，就是假設，比較，修正。當前，迭代圖像重建算法有兩類，即代數(shù)迭代重建算法與統(tǒng)計迭代重建算法。代數(shù)迭代重建算法是以代數(shù)方程理論為基礎，主要有一般ART算法和同時代數(shù)重建算法。而統(tǒng)計迭代重建算法是基于各種統(tǒng)計準則，主要有最小均方誤差、最大似然估計等重建算法。3.3.1代數(shù)迭代重建算法(1)ART算法ART算法有時也稱為Kaczmarz算法，它的主要思路是讓當前所有估算的圖像在每一次更新中滿足一個方程，在迭代修正過程中，每次只考慮一個投影單元的投影值。原理如下圖3-7所示。比較并計算修正量比較并計算修正量理論投影值實際投影值斷層圖像估計值迭代循環(huán)環(huán)環(huán)環(huán)修改計算3-7迭代算法計算過程示意圖最常見的ART算法是基于交替投影法進行迭代修正的，它的圖像更新式為： (3-19)為松弛參數(shù)。(2)同時代數(shù)重建算法由于ART算法對圖像值進行修正時只依賴一條投影帶上數(shù)據(jù)，因此人們又提出了同時代數(shù)重建算法（SART算法），它是在校正像素單元的圖像值之前，計算出像素單元上所有投影估計值與實際值的差別，并求出來，再利用平均值對圖像進行修正。SART算法的迭代公式如下： (3-20)這樣，經過各條投影帶上的平均值，能夠減小誤差，避免對重建結果帶來過大影響，同時它又抑制圖像重建過程中的噪聲。3.3.2影響代數(shù)迭代重建算法的因素ART重建算法之因此不能廣泛應用，是因為它是一個重復迭代與修正的過程，計算量大。經過長期探索與研究，得到了影響ART算法的很多因素。如下：(1)基函數(shù)的選擇；(2)松弛參數(shù)的選擇；(3)迭代次數(shù)的最優(yōu)設計；(4)數(shù)據(jù)訪問方式?；瘮?shù)的選擇直接影響重建算法的計算量。如果基函數(shù)選擇過于復雜，重建的圖像更加逼近真實的，可是計算量太龐大，如果選擇很簡單，雖然計算量小，但重建圖像質量不高。ART重建算法就是解線性方程并伴隨松弛參數(shù)的迭代過程，這樣選擇合理的松弛參數(shù)會很大程度地影響運算速度?，F(xiàn)如今，松弛參數(shù)的選擇的基本標準是。只有在這個區(qū)間迭代重建算法才能按照標準并保證它的收斂性[6]。雖然ART迭代重建是不斷重建的過程，可是它的重建質量不與迭代次數(shù)成正比，而是在達到某個迭代次數(shù)后，質量就開始下降。因此選擇最優(yōu)的迭代次數(shù)是必要的，不但節(jié)省時間又有高質量的重建圖像。下面我們用一個實例來說明這種情況。我們選用一個圓形工件，該圓形工件內部含有碳、氣泡、鉛等雜質，工件的一些雜質參數(shù)[7]如下表3.1所示。名稱半徑(cm)圓心坐標(x,y)吸收系數(shù)(1/cm)鐵5.0(0.0,0.0)0.435碳2.0(0.0,2.5)0.115空氣0.5(0.0,-1.0)0.000鉛1.0(0.0,-3.0)0.800表3.1工件的雜質參數(shù)表根據(jù)參數(shù)，我們到模型的二維灰度分布圖及y軸上一維分布圖如下圖3.8所示。3-8模型的二維灰度分布圖及y軸上一維分布圖選擇松弛因子為0.1,抽樣間距為0.1cm，重建的像素為00的重建圖像，得到二維灰度分布與y軸上一維分布圖如下圖3-9至3-13所示。3-8迭代次數(shù)為1次3-9迭代次數(shù)為3次3-10迭代次數(shù)為5次3-11迭代次數(shù)為10次3-12迭代次數(shù)為15次經過上述圖片，我們能夠看出第1到第3次迭代圖像的質量有明顯的改進，第5到第10次的迭代圖像質量緩慢變好，從第10次到第15次迭代圖像質量明顯下降。因此，并不是一味的增加迭代次數(shù)就能夠改進重建圖像質量。3.3.3ART重建算法與SARTART算法每次迭代時只有一條投影值，如果在這條射線在投影時引入了噪聲，這樣對重建圖像也會帶來誤差。而SART重建算法的每個像素校正值是所有射線的誤差累計和，這樣它能有效地抑制測量數(shù)據(jù)中的噪聲。SART算法的迭代速度比ART算法快很多，大約差一個數(shù)量級。ART在選擇松弛因子方面范圍較大，在(0,2)之間，而SART算法在(0,0.001)之間，選擇性較窄。盡管SART算法抗噪性較好，可是由于在選擇松弛因子方面比較小，同樣使得重建圖像很模糊，如下圖3-13和3-14所示，在高斯隨機噪聲為情況下，迭代次數(shù)為5次的條件下，ART與SART算法的重建圖像。3-13松弛因子為0.1迭代次數(shù)為5次的ART重建圖像3-14松弛因子為0.0001迭代次數(shù)為5次的SIRT重建圖像從上圖看出，雖然SART算法重建圖像比較光滑，可是圖像比較模糊。3.3.4統(tǒng)計迭代重建算法(1)OS-EM算法在OS-EM（分成有序的子集來求期望值的極大值）算法中，數(shù)據(jù)被分成不同的子集，按照給定的次序走訪不同的子集，每訪問一個子集圖，像就更新一次[8]。子集的分割有很大的自由度，能夠根據(jù)需要進行調整，不必按照投影角度來劃分。一般來說，若使用N個子集，收斂速度大約提升N倍?？墒歉鶕?jù)人們的經驗，一般讓子集分割達到收斂速度提升10倍左右，在這樣條件下噪聲破壞不明顯。第四章濾波反投影重建算法濾波反投影算法（FBP）是CT圖像重建中的一種經典算法，因為它的計算效率高，需求資源少，因此被廣泛地應用在各個領域上，特別是醫(yī)學領域。反投影重建算法是直接對投影進行重建并得出圖像，這種方法重建得出的圖像質量不高，容易失真。與反投影重建算法不同的是，濾波反投影算法是先對投影數(shù)據(jù)進行濾波，再對數(shù)據(jù)進行反投影重建，這樣得出的圖像更加清晰準確。根據(jù)成像系統(tǒng)采集數(shù)據(jù)方式來分，又可分為平行束濾波反投影重建算法和扇形束濾波反投影算法[9]。FBP重建算法有兩種計算方法，一種與中心切片定理有密切關系，叫做卷積反投影重建算法。另一種是Radon反變換，我們常見第一種，下面主要對第一種做詳細介紹。4.1平行束濾波反投影重建算法4.1.1卷積反投影重建算法卷積反投影算法是以中心切片定理為基礎產生的。由中心切片定理知，重建圖像的二維傅里葉變換可由在不同視角下的投影的，F(xiàn)BP重建算法坐標示意圖如下圖4-1所示。4-1FBP算法坐標系需要重建圖像為： (4-1)如圖所示，，，，能夠得到 (4-2)經過上式變換可得到， (4-3)式中，,,。事實上，圖像重建的過程能夠看成是由一系列坐標變換得到的。由上列推倒知，F(xiàn)BP算法是在一定視角下投影，然后進行濾波投影，再做反投影，把這些反投影值累加就能夠得到重建圖像。4.1.2濾波函數(shù)的選擇濾波函數(shù)的選擇在FBP算法中扮演一個重要的角色。理想的濾波函數(shù)能夠使重建圖像更加準確清晰。上面式子中的濾波函數(shù)是理想化的，不是平方可積的，因此濾波函數(shù)無法實現(xiàn)。合理地選擇濾波函數(shù)是必要的，常見的主要有R-L和S-L兩種濾波函數(shù)。(1)R-L濾波函數(shù)系統(tǒng)函數(shù)表示為： (4-4)其中，代表采樣間隔，為最高不失真頻率。系統(tǒng)沖激函數(shù)表為 (4-5)我們在對圖像濾波時，常常見離散化的函數(shù)，這樣將帶入上述式子，這樣離散化的沖激函數(shù)表示為： (4-6)函數(shù)的連續(xù)空域特性如下圖4-2所示。4-2沖激函數(shù)空域特性經過上式，我們能夠看到它形式簡單，實用性強，用它來重建的圖像輪廓上清楚?？墒且灿腥秉c，有吉布斯效應，表現(xiàn)為明顯的振蕩響應。另外，重建的圖像容易受到噪聲的影響。(2)S-L濾波函數(shù)此函數(shù)緩解了R-L濾波函數(shù)的不足，選擇了不同的窗函數(shù),所對應的濾波函數(shù)也不相同。該濾波函數(shù)所對應的系統(tǒng)函數(shù)為： (4-7)這里的以及的含義與R-L濾波器中的一樣。它的沖激響應函數(shù)為： (4-8)同上述R-L濾波函數(shù)一樣，可得到沖激函數(shù)離散化形式為： (4-9)函數(shù)的空域特性如下圖4-3所示。4-3函數(shù)空域特性利用S-L濾波函數(shù)來重建圖像的優(yōu)點是，振蕩響應小，而且在投影數(shù)據(jù)有噪聲的條件下，它的重建質量也較R-L好，可是在低頻段上，R-L濾波函數(shù)的重建質量相對較高[10]。4.2扇形束濾波反投影算法由于扇形束CT系統(tǒng)探測器陣列方式有圓弧形和直線型，這樣扇形射線就分為等角型和等間距型。相應的重建算法就有等夾角型和等間距型扇形束濾波反投影算法。如下圖所示4-4和4-5所示。4-4等夾角型FBP算法4-5等間距型FBP算法下面將對兩種方法分別討論。4.2.1等夾角扇形束濾波反投影算法等夾角扇形束的參數(shù)關系圖如下圖4-6所示。4-6等夾角扇形束參數(shù)關系圖其中，表示為扇形區(qū)域內部物體，為射線源，弧線為探測器所在位置，為射線源到中心的距離，為中心射線，扇形的位置由和夾角決定，任一條射線由決定，由于坐標系固定，因此射線位置就取決于。設投影函數(shù)為，經過投影函數(shù)重建圖像。扇形束的FBP重建算法可由平行束的FBP重建算法推導出來。首先，將公式變換得到： (4-10)其中，，,為物體斷層內任一點的極坐標。經過上式能夠看出，扇形束的FBP重建算法是對平行束FBP重建算法的加權與修正。4.2.2等間距扇形濾波反投影算法等間距扇形束參數(shù)關系圖如圖4-7所示。4-7等間距扇形束參數(shù)關系圖同樣的，是扇形區(qū)域內部的物體，線段為探測器所在位置，為任意射線，為射線源與中心距離，由于坐標系固定，這樣由確定。同樣是虛擬探測器。設投影函數(shù)為，現(xiàn)利用關系求重建圖像。同等夾角扇束FBP重建表示式推導過程一樣，由于平行束的數(shù)據(jù)與扇形束的數(shù)據(jù)不同，只需將一些變量修改，經過變換，最后得到的等間距扇束FBP重建算法的表示式為： (4-11)其中，是等效投影，，為斷層內任一點極坐標。4.3濾波反投影重建算法的軟件實現(xiàn)這里，我們選擇MATLAB與VisualC++相結合的方式對圖像進行重建。我們知道，VisualC++軟件更注重類的創(chuàng)立，而且經過在每個類中編寫相應類的程序能更加深入地了解圖像重建的過程，然而對于MATLAB軟件來說，因為它本身具有相應的函數(shù)，因此利用它能簡潔地分析參數(shù)對重建圖像的影響。首先我們選用經典的shepp-logan頭部模型來研究，該模型是由11個不同密度的橢圓疊加而成的，這種模型被用來進行圖像的仿真與評估算法的性能。在VisualC++軟件中，根據(jù)面向對象的思想，在photostar平臺中創(chuàng)立CCTEmulate類來實現(xiàn)仿真。實驗中，讀入圖像為灰度圖像即可，在CCTEmulate類中設置模型函數(shù)的菜單消息，要確定像素在哪個橢圓中。Shepp-logan函數(shù)的思路是：檢測圖像的每個像素是否依次在某一個橢圓中，若在橢圓旋轉角度外，置為背景色，若在某個橢圓內，置為其灰度值。這些操作都在歸一化坐標系中完成。在相應類中，編寫程序如下。BOOLCCTEmulate::SheepLogan(){ ASSERT(m_pimageproject!=NULL); intnWidth=pImageObject->GetWidth(); intnHeight=pImageObject->GetHeight(); intnNumBits=pImageObject->GetNumBits(); intnWidthBytes; char*pBuffer=(char*)pImageObject->GetDIBPointer(&nWidthBytes); if(pBuffer==NULL)return(FALSE); BITMAPFILEHEADER*pBFH; BITMAPINFOHEADER*pBIH; RGBQUAD*pRGBPalette; unsignedchar*pBits; intnNumColors=pImageObject->GetNumColors(); pBFH=(BITMAPFILEHEADER*)pBuffer; pBIH=(BITMAPINFOHEADER*)&pBuffer[sizeof(BITMAPFILEHEADER)]; pRGBPalette=(RGBQUAD*)&pBuffer[sizeof(BITMAPFILEHEADER) +sizeof(BITMAPINFOHEADER)]; pBits=(unsignedchar*)&pBuffer[sizeof(BITMAPFILEHEADER) +sizeof(BITMAPINFOHEADER)+nNumColors*sizeof(RGBQUAD)]; for(inti=0;i<nWidth;i++) {for(intj=0;j<nHeight;j++) {intk=0; switch(k) { case0://橢圓a if(EllipseCaculate(0,0,0.92,0.69,0,1,i,j)>1) { pBits[(nWidth-1-i)*nWidthBytes+j]=0*255;break; } elsepBits[(nWidth-1-i)*nWidthBytes+j]=1*255;case1://橢圓b if(EllipseCaculate(0,-0.0184,0.874,0.6624,0,1,i,j)>1)break;elsepBits[(nWidth-1-i)*nWidthBytes+j]=0.5*255; case2://橢圓cif(EllipseCaculate(0.22,0,0.31,0.11,0.3090,0.9511,i,j)<=1) { pBits[(nWidth-1-i)*nWidthBytes+j]=int(0.3*255); break;以下橢圓本分省略。}在產生投影數(shù)據(jù)時，要確定合理的投影數(shù)，即在歸一化坐標系中選擇理想的角度進行投影，這里我們將圖像的指定直線像素置為255，每10度畫線一次，根據(jù)要求編寫程序。 intnProjectionHalfLength=(nWidth>nHeight)?nWidth:nHeight; nProjectionHalfLength=(int)((nProjectionHalfLength+1)/2.0);intnProjectionLength=2*nProjectionHalfLength-1;for(inti=0;i<180;i++)//進行投影，每一度一次 { floatr1=i/180.0*3.1415;//把角度轉換為弧度 floatr2=(i+180)/180.0*3.1415; for(intk=0;k<nProjectionLength;k++) {//法線長度 floatL=float(k-nProjectionHalfLength+1);//原圖像的列,從0開始計算 intn=0;while(k<nProjectionHalfLength-1&&n<nWidth)//r2 {floatx=n-X; floaty=(L-x*cos(r2))/sin(r2); intm=int(Y-y);n++; if(m<0||m>=nHeight)continue; pProjDa[i*nProjectionLength+k]=pProjDa[i*nProjectionLength+k]+pBits[(nHeight-1-m)*nWidthBytes+n];} n=0;//對每一條直線x軸都要掃描一次 while(k>=nProjectionHalfLength-1&&n<nWidth)//r1 { floatx=n-X; floaty=(L-x*cos(r1))/sin(r1); intm=int(Y-y);開始計算 n++; if(m<0||m>=nHeight) continue;pProjDa[i*nProjectionLength+k]=pProjDa[i*nProjectionLength+k]+pBits[(nHeight-1-m)*nWidthBytes+n]; }選擇濾波反投影重建圖像時，需要解決濾波器的采樣頻率，從空間坐標系與歸一化坐標系的關系看，濾波器采樣頻率在圖像采樣頻率的2倍以上即可。我們接下來進行構造濾波器，編寫相應的程序。根據(jù)以上基本思想，我們編寫程序圖，進入界面如下圖4-8所示。4-8編寫程序運行圖運行之后如下圖4-9所示，4-9生成TEST文本圖選擇仿真按鈕，生成shepp-logan模型圖片如下4-10所示，4-10shepp-logan模型圖4-11CT掃描過程圖對shepp-logan模型圖片進行重建仿真，得到過程圖片如上圖4-11所示。重建仿真之后，選擇濾波反投影成像之后，我們能夠看到有兩種濾波函數(shù)，一種是R-L濾波器，另一種是S-L濾波器。如下圖4-12所示。4-12兩種濾波函數(shù)界面選擇R-L濾波器進行仿真，得到圖像如下圖4-13所示。4-13R-L濾波重建圖像4-14S-L濾波重建圖像選擇S-L濾波器進行仿真，得到圖像如上圖4-14所示。經過濾波反投影重建圖像之后，我們能夠看到，圖像能大致反應原來圖像的特性，可是分辨率比較低，原始圖像中的下三個橢圓幾乎無法形成，圖像的上下邊緣有大面積的白色偽跡。相比于R-L濾波器，S-L濾波器重建效果要好些，白色雜點的狀況相對較少[11]。其實我們能夠用MATLAB仿真軟件進行，因為MATLAB里面含有自帶的函數(shù)，只需調整相應參數(shù)即可。下面用MATLAB來實現(xiàn)圖像重建，并分析投影數(shù)對重建圖像的影響。用MATLAB中圖像處理工具箱中的phantom生成shepp-logan模型。程序：P=phantom(256);Imshow(P)用MATLAB中的radon函數(shù)獲得shepp-logan模型的投影數(shù)據(jù)。程序：thetal1=0:10:170;[R1,xp]=radon(P,thetal);thetal2=0:5:175;[R2,xp]=radon(P,theta2);thetal3=0:2: