第一章 緒論

1.1 引言

計(jì)算機(jī)斷層攝像技術(shù)CT(Computerized Tomography)是20世紀(jì)醫(yī)學(xué)的重大成果之一，該成果將計(jì)算機(jī)應(yīng)用于醫(yī)學(xué)領(lǐng)域，使醫(yī)學(xué)射線學(xué)發(fā)生了革命性的變化。自從CT問世以來，計(jì)算機(jī)科學(xué)與生物醫(yī)學(xué)工程相結(jié)合，形成了計(jì)算機(jī)醫(yī)學(xué)圖像研究的新領(lǐng)域，并為生命科學(xué)的研究提供了新的方法，成為近年來世界科技界最活躍、最富有生機(jī)和成就的領(lǐng)域之一。

醫(yī)學(xué)影像學(xué)將數(shù)字圖像處理技術(shù)和計(jì)算機(jī)圖形學(xué)技術(shù)廣泛的應(yīng)用于生物醫(yī)學(xué)領(lǐng)域中，通過把人體的內(nèi)部結(jié)構(gòu)以圖像或圖形的方式顯示出來，提高了醫(yī)療診斷的可靠性，使治療能夠準(zhǔn)確和徹底。

1.2 醫(yī)用CT的簡介

CT是計(jì)算機(jī)X射線斷層造影術(shù)(Computerized Tomography)的簡寫。CT的發(fā)明是20世紀(jì)后期最重大的科技成果之一，由Hounsfield于1969年設(shè)計(jì)成功，1972年公諸于世。

CT利用人體各種組織(包括正常和異常組織)對X射線的吸收不等這一特性，將人體某一選定層面分成許多立方體小塊(這些立方體小塊稱為體素)，X射線通過人體測得每一體素的密度或灰度，即為CT圖像上的基本單位，稱為像素。它們排列成行列方陣，形成圖像矩陣。當(dāng)X射線球管從一方向發(fā)出X射束穿過選定層面時(shí)，沿該方向排列的各體素均在一定程度上吸收一部分X射線，使X射線衰減。當(dāng)該X射線束穿透組織層面(包括許多體素)為對面探測器接收時(shí)X射線量已衰減很多，為該方向所有體素X射線衰減值的總和。然后X射線球管轉(zhuǎn)動(dòng)一定角度，再沿另一方向發(fā)出X射線束，則在其對面的探測器可測得沿第2次照射方向所有體素X射線衰減值的總和;以同樣方法反復(fù)多次在不同方向?qū)M織的選定層面進(jìn)行X射線掃描，即可得到若干個(gè)X射線衰減值總和。在上述過程中，每掃描一次，即可得一方程。該方程中X射線衰減總量為已知值，而形成該總量的各體素X射線衰減值是未知值。經(jīng)過若干次掃描，即可得一聯(lián)立方程組，經(jīng)過計(jì)算機(jī)運(yùn)算可解出這一聯(lián)立方程組，而求出每一體素的X射衰減值，再經(jīng)模/數(shù)轉(zhuǎn)換，使各體素不同的衰減值形成相應(yīng)各像素的不同灰度，各像素所形成的矩陣圖像即為該層面不同密度組織的灰度圖像。

螺旋CT檢查包括兩方面的基本內(nèi)容：一是X射管及探測器連續(xù)360°旋轉(zhuǎn);二是患者同時(shí)隨檢查床勻速推進(jìn)，如圖1.2所示。在掃描時(shí)間內(nèi)，X射線焦點(diǎn)對病人作螺旋式運(yùn)動(dòng)，并同時(shí)收集這一范圍的全部掃描數(shù)據(jù)，用線性內(nèi)插法重建圖像。

如圖1.3所示，醫(yī)用X-CT機(jī)的系統(tǒng)結(jié)構(gòu)主要由六大部分組成，其各部分的作用如下：

(1)X射線源：產(chǎn)生用來檢測被測物的X射線，X射線源包括X射線球管源(能量在450kV以下)和直線加速器(能量在2MeV以上)。射線源的能量，決定了穿透能力。

(2)探測系統(tǒng)：包括準(zhǔn)直器、傳感器、信號(hào)處理和信號(hào)傳輸?shù)炔糠?，是獲取信號(hào)的關(guān)鍵部分，也是決定CT性能的關(guān)鍵部分之一。穿過被測物的X射線首先通過準(zhǔn)直器準(zhǔn)直并離散化，傳感器先將射線轉(zhuǎn)換成電信號(hào)，信號(hào)處理電路再將不規(guī)則的信號(hào)轉(zhuǎn)換成標(biāo)準(zhǔn)的信號(hào)傳輸?shù)接?jì)算機(jī)接口。

(3)計(jì)算機(jī)采集系統(tǒng)：主要由特殊的專用的多信道數(shù)據(jù)采集接口電路和計(jì)算機(jī)軟硬件組成。完成數(shù)據(jù)采集、轉(zhuǎn)換、校正、處理等。將采集的數(shù)據(jù)處理成標(biāo)準(zhǔn)的文件格式，供圖像重建、處理使用。

(4)機(jī)械掃描系統(tǒng)：作為各部分的載體并提供CT掃描所需的多個(gè)自由度的高精度運(yùn)動(dòng)。

(5)自動(dòng)控制系統(tǒng)：包括檢測、驅(qū)動(dòng)、控制器(計(jì)算機(jī))，完成掃描運(yùn)動(dòng)控制、系統(tǒng)邏輯和程控、狀態(tài)監(jiān)控和安全保護(hù)，協(xié)調(diào)整機(jī)工作，并完成系統(tǒng)自檢與數(shù)據(jù)診斷。

(6)圖像處理系統(tǒng)：包括圖像處理計(jì)算機(jī)硬件和軟件，如用于圖像重建與處理的高速計(jì)算機(jī)、大屏幕圖像顯示器、大容量數(shù)據(jù)存貯器、圖像拷貝輸出設(shè)備(打印機(jī))、系統(tǒng)軟件及專用軟件。完成數(shù)據(jù)校正、圖像重建、處理、分析、測量、圖像輸出、存貯、顯示等。

我們所研究的CT圖像重建部分處于圖像處理系統(tǒng)中，是整個(gè)系統(tǒng)的瓶頸所在，也是決定系統(tǒng)整個(gè)過程所消耗時(shí)間的關(guān)鍵部分。

1.3 CT圖像重建技術(shù)概述

1.3.1 CT圖像重建的簡介

我們試圖重建的物體可被看作是某種函數(shù)的二維分布。對于CT，該函數(shù)代表物體線性衰減系數(shù)。關(guān)于斷層重建問題的描述，我們可以假設(shè)采集了一組測量結(jié)果，每個(gè)測量結(jié)果代表沿著特定的射線路徑，物體衰減系數(shù)的累加或線積分。這些測量結(jié)果是在不同角度和到旋轉(zhuǎn)中心的不同距離上獲取的。為避免數(shù)據(jù)采樣的冗余，我們假設(shè)測量按以下次序進(jìn)行。首先沿著彼此平行且等間距的路徑進(jìn)行一組測量。這些測量結(jié)果構(gòu)成一次“觀測”或一組“投影”。在略微改變的角度下重復(fù)同樣的測量。持續(xù)該過程直到覆蓋整個(gè)360°(理論上僅有180°平行投影是必要的)。在整個(gè)過程中，相鄰兩次觀測之間的角度增量保持不變，并且被掃描物體在同一位置固定不動(dòng)。CT重建的問題就是，我們?nèi)绾位谶@些測量結(jié)果來估計(jì)被掃描物體的衰減系數(shù)分布。

CT圖像重建問題是一個(gè)有趣而復(fù)雜的課題。它的公示表達(dá)可以追溯到1917年，當(dāng)時(shí)Radon(雷登)[2]首先找到了從函數(shù)線積分重建該函數(shù)的求解方法。隨著20世紀(jì)70年代后期和80年代早期臨床實(shí)用CT掃描機(jī)的發(fā)展，該領(lǐng)域的研究活動(dòng)有了極大的發(fā)展。大量研究論文、會(huì)議論文匯編、書籍章節(jié)，甚至教科書都關(guān)注這個(gè)課題[3,4]。提出了許多技術(shù)，它們在計(jì)算復(fù)雜性、空間分辨率、時(shí)間分辨率、噪聲、臨床治療方案、靈活性以及偽像各方面具有不同的折中平衡。

1.3.2 Radon(雷登)變換

CT的基本思想源于1917年奧地利數(shù)學(xué)家Radon提出的Radon變換。

Radon變換的內(nèi)容可以表述為：若已知某函數(shù)，

如圖1.4所示，其沿直線S的線積分為：

(1.1)

則

(1.2)

式(1.1)為Radon變換，實(shí)際上就是物體的投影，式(1.2)為Radon反變換，即根據(jù)投影數(shù)據(jù) 重建函數(shù) 。

1.3.3 傅里葉切片定理

傅里葉切片定理的含義是：平行投影的一維傅里葉變換等同于原始物體的二維傅里葉變換的一個(gè)切片。即是指出線性衰減系數(shù)函數(shù)f(x,y)在某一方向上的投影函數(shù)gθ(R)的一維傅立葉變換函數(shù)Gθ(ρ)是f(x,y)的二維傅立葉變換函數(shù)F(u,v)或F(ρ,θ)(極坐標(biāo)形式)在(ρ,θ)平面上沿同一方向且過原點(diǎn)的直線上的值，如圖1.5所示。

為此，我們在不同的角度下取得足夠多的投影函數(shù)數(shù)據(jù)，并作它們的傅立葉變換，那么變換后的數(shù)據(jù)就將充滿整個(gè)(u,v)平面。一旦頻域函數(shù)F(u,v)或F(ρ,θ)的全部值得到后，將其作一次傅立葉反變換，就得到原始的衰減系數(shù)函數(shù)f(x,y),即

(1.3)

令u=ρcosθ，v=ρsinθ，則式(1)可進(jìn)一步變形為

(1.4)

式中 ，表示對投影函數(shù) 的傅里葉變換函數(shù)進(jìn)行濾波變換，其中 為濾波函數(shù)。

由傅立葉變換性質(zhì)可知，頻域中的濾波運(yùn)算可等效地在空域中用卷積運(yùn)算來完成，因此由(2)可得到

(1.5)

式中h(R)為濾波函數(shù) 的空域形式， ，因而這種方法也稱為卷積反投影方法。

利用中心切片定理[5]及二維FFT反變換法重建圖像，由于勿須反投影運(yùn)算,因而速度快,但圖像重建過程中,需要內(nèi)插運(yùn)算，因而重建圖像精度相對較低。

首先求出各投影數(shù)據(jù)的一維傅里葉變換，在不同的投影角度下所到的一維變換函數(shù)可構(gòu)成完整的二維傅里葉變換函數(shù)，將此二維函數(shù)作一次反傅里葉變換，就得到重建圖像。為了在二維逆變換中采用快速傅里葉變換算法，通常在逆變換前要將極坐標(biāo)轉(zhuǎn)化為直角坐標(biāo)的形式。

傅里葉變換法重建法的特點(diǎn)是變換速度快，但精度不如濾波反投影法。算法的關(guān)鍵是將弧形的的極坐標(biāo)數(shù)據(jù)轉(zhuǎn)換成直角坐標(biāo)數(shù)據(jù)時(shí)，由于在邊緣區(qū)高頻數(shù)據(jù)減少，因而造成誤差，但傅里葉變換重建法重建速度比濾波反投影可提高2-3倍一在弧形極坐標(biāo)數(shù)據(jù)向直角坐標(biāo)系轉(zhuǎn)化時(shí)，最簡單的是最鄰近內(nèi)插法，當(dāng)然這種方法精度最低，雙線性內(nèi)插重建圖像精度好于最鄰近內(nèi)插法，而且計(jì)算又不復(fù)雜。

解決的方法是擴(kuò)大計(jì)算區(qū)域，通過外延數(shù)據(jù)附加上一些格外的點(diǎn)，即計(jì)算更多的像素點(diǎn)以減小邊緣的誤差。如重建圖像為M×M，則可計(jì)算3M×3M區(qū)域內(nèi)的FFT變換，當(dāng)然這是以增加了計(jì)算量為代價(jià)的。傅里葉變換重建圖像算法在內(nèi)插網(wǎng)格點(diǎn)上進(jìn)行一些適當(dāng)?shù)倪x擇。如使徑向點(diǎn)取在直角坐標(biāo)網(wǎng)格的線上，這樣只需一次內(nèi)插，而重建圖像精度有了較大的改進(jìn)。

1.3.4 CT圖像重建的幾種算法

在實(shí)際重建當(dāng)中所存在的問題是，雖然Radon給出了一個(gè)數(shù)學(xué)公式，但是我們需要一個(gè)有效的算法來解決它，圖像重建的算法有很多，大致分為三類：精確算法、近似算法和迭代算法。近似算法中，以濾波反投影算法(Filter back projection，F(xiàn)BP)最具代表性，應(yīng)用最為廣泛。選代算法中，代數(shù)重建算法(Algebraic reconstruction technique，ART)是提出最早并最為人們熟悉的算法。迭代型算法(如代數(shù)重建算法等)具有許多優(yōu)點(diǎn)，但由于計(jì)算量大、重建時(shí)間長.在很長一段時(shí)間內(nèi)限制了其在醫(yī)學(xué)和工業(yè)CT領(lǐng)域的應(yīng)用。提高迭代型算法的計(jì)算速度一直是人們關(guān)注的問題。近年來人們提出了不少提高迭代算計(jì)算速度的方法，加上近年來計(jì)算機(jī)計(jì)算速度的迅速提高，迭代算法重新受到人們青睞。此外，由于應(yīng)用的需要，局部重建算法(Local Reconstruction Algorithm， LocalRA)也在近十年中有了較大的發(fā)展。在傳統(tǒng)全局CT算法中，即使重建物體斷面中一個(gè)小區(qū)域的圖像，也得圍繞整個(gè)斷面采集投影數(shù)據(jù)。而局部重建算法，僅需圍繞感興趣區(qū)域及其鄰域采集投影數(shù)據(jù)，即可重建感興趣區(qū)域的圖像。局部重建算法可減少數(shù)據(jù)采集時(shí)間和重建時(shí)間，降低人體(或生物體)的放射攝入量。[!--empirenews.page--]

1.3.5國內(nèi)外研究現(xiàn)狀

同類課題所研究的技術(shù)基本上被國外所壟斷，國內(nèi)尚未有人提出，國內(nèi)現(xiàn)在所使用的技術(shù)是利用PC機(jī)上軟件來實(shí)現(xiàn)圖像的重構(gòu)，所需時(shí)間較長，如果用FPGA來實(shí)現(xiàn)的話，速度可以提高數(shù)十乃至上百倍。

1.4研究背景及意義

在當(dāng)今社會(huì)大力發(fā)展醫(yī)療衛(wèi)生條件的背景下，許多醫(yī)院迫切需要先進(jìn)的CT來為患者診斷病情，現(xiàn)在的CT技術(shù)被國外所壟斷，設(shè)備也都在200萬以上，只有極少數(shù)醫(yī)院有能力配備，所以急需研發(fā)具有自主知識(shí)產(chǎn)權(quán)的產(chǎn)品，把價(jià)格控制在50萬以內(nèi)。CT的關(guān)鍵技術(shù)之一是快速斷層圖像重建技術(shù)，本課題的立足點(diǎn)就在于利用FPGA的高度并行性，實(shí)現(xiàn)CT斷層圖像重建算法，滿足實(shí)際產(chǎn)品速度要求，為實(shí)現(xiàn)CT國產(chǎn)化準(zhǔn)備，推動(dòng)社會(huì)醫(yī)療衛(wèi)生條件的發(fā)展。

第二章 濾波反投影算法

2.1 濾波反投影算法介紹

盡管傅里葉切片定理提供了斷層成像重建的一個(gè)直接方案，在真正實(shí)現(xiàn)過程中，它提出了一些難題。首先，傅里葉空間中產(chǎn)生的采樣模式不是笛卡兒坐標(biāo)的。傅里葉切片定理說明一次投影的傅里葉變換是二維傅里葉空間中通過原點(diǎn)的一條直線。結(jié)果，不同投影采樣落到極坐標(biāo)柵格上。為了執(zhí)行二維傅里葉變換，這些采樣不得不被插值或重新柵格化到一個(gè)笛卡兒坐標(biāo)中。二維頻率域中的插值不像真實(shí)空間中的插值一樣直接。在真實(shí)空間里，一個(gè)插值誤差局限于像素所在的小區(qū)域。然而，對于頻域插值，這個(gè)特性不再有效，因?yàn)槎S傅里葉空間中每個(gè)采樣表示某一個(gè)空間頻域(在水平和垂直方向上)。于是，在傅里葉空間中一個(gè)單獨(dú)采樣點(diǎn)上產(chǎn)生的誤差會(huì)影響整個(gè)圖像(經(jīng)過傅里葉反變換后)的外貌。為闡明傅里葉域插值的敏感性，進(jìn)行下面的簡單實(shí)驗(yàn)。掃描一個(gè)肩部模體，并在512×512矩陣中重建，矩陣用f(x,y)表示，其中x=0，1，…512，y= 0，1，…，512。下一步，執(zhí)行圖像的二維離散傅里葉變換，得到一個(gè)函數(shù)F(u,v)，其中u=0，1，…，511，v=0，1，…，511。注意F(u,v)是一個(gè)512×512復(fù)數(shù)矩陣。在該矩陣中，F(xiàn)(00)代表圖像的直流成分。如果簡單地進(jìn)行函數(shù)F(u,v)的離散傅里葉反變換，將得到原始圖像f(x,y)。注意函數(shù)F(u,v)是我們試圖采用平行投影進(jìn)行估計(jì)的量值(傅里葉切片定理)。

直接傅里葉域重建的另一缺點(diǎn)是進(jìn)行目標(biāo)重建的困難性。目標(biāo)重建是在CT中常用的技術(shù)，用來檢查物體中一個(gè)小區(qū)域的精密細(xì)節(jié)。如果能以某種方式把重建“聚焦”在感興趣區(qū)，物體的細(xì)節(jié)就可以更好地顯現(xiàn)。采用直接傅里葉重建方法，需要用大量的0填充F(u,v)，以進(jìn)行必要的頻率域插值。傅里葉反變換的大小和目標(biāo)ROI的尺寸成反比。對非常小的ROI，矩陣尺寸龐大以至無法管理。盡管其他技術(shù)可以用來克服其中一些困難，這些技術(shù)的實(shí)現(xiàn)仍不直截了當(dāng)。因此，必須研究傅里葉切片定理的替代實(shí)現(xiàn)方法。濾波反投影算法是目前得到廣泛應(yīng)用的基于變換法的圖像重建算法，它具有重建速度快、空間和密度分辨率高等優(yōu)點(diǎn)，缺點(diǎn)是對投影數(shù)據(jù)的完備性要求較高[7]，從數(shù)學(xué)上講，只有獲得被檢試件所有的Radon變換數(shù)據(jù)(完全投影數(shù)據(jù))后才能精確重建其切片圖像。

2.2 濾波反投影算法公式的推導(dǎo)

我們從傅里葉變換和傅里葉反變換是共扼算子這一眾所周知的事實(shí)開始。圖像函數(shù)f(x,y)可以通過傅里葉反變換從它的傅里葉變換F(u,v)中恢復(fù)，

(2.1)

與推導(dǎo)傅里葉切片定理時(shí)進(jìn)行的坐標(biāo)變換類似，我們從笛卡兒直角坐標(biāo)(u,v)轉(zhuǎn)換到極坐標(biāo)

。

坐標(biāo)轉(zhuǎn)換的目的是以更自然的數(shù)據(jù)采集形式表達(dá)數(shù)值F(u,v)。坐標(biāo)轉(zhuǎn)換如下：

，

， (2.2)

和

(2.3)

將等式(2.2)和(2.3)帶入到(3.1)，得到

(2.4)

利用公式

中描述的傅里葉切片定理，我們用代替，建立如下關(guān)系：

(2.5)

.

對于平行采樣幾何束，在投影采樣中存在一個(gè)微妙的對稱性：

(2.6)

通過研究一組相差180°的平行束的采樣幾何，這個(gè)特性可以很容易理解。兩組投影正好代表同一組射線路徑。基于傅里葉變換的特性，對于相應(yīng)的傅里葉變換對來說，存在一個(gè)簡單關(guān)系：

(2.7)

將等式(2.7)代入等式(2.5)，我們得到下面等式：

(2.8)

通過在旋轉(zhuǎn)坐標(biāo)系(s,t)中表達(dá)上面的等式，并利用等式：

.

中指出的關(guān)系，我們得到下面等式：

. (2.9)

這里，是在角度投影的傅里葉變換。內(nèi)部積分是數(shù)值的傅里葉反變換。在空間域，它代表一個(gè)經(jīng)頻域響應(yīng)為的函數(shù)濾波后的投影。我們稱之為“濾波投影”。

如果用標(biāo)記等式(2.9)的內(nèi)部積分所代表的角上的濾波投影：

. (2.10)

等式(2.8)可以下面形式重寫

(2.11)

變量是從點(diǎn)(x,y)到一條通過坐標(biāo)系原點(diǎn)，并與x軸成角的直線的距離。等式(2.11)說明，重建圖像f(x,y)在位置(x,y)，是通過該點(diǎn)的所有濾波投影采樣的累加。另外，我們還可以選擇關(guān)注一個(gè)特定濾波投影采樣，研究它對重建圖像的貢獻(xiàn)。因?yàn)榇砼c產(chǎn)生投影采樣的射線路徑重疊的一條直線，的強(qiáng)度沿著直線均勻地加到重建圖像。結(jié)果，濾波投影采樣的值沿著整個(gè)直線路徑被“涂抹”或“疊加”。

我們還可以給出濾波反投影方法的一個(gè)直觀解釋?；诟道锶~切片定理，物體的二維傅里葉變換是通過將許多一維傅里葉變換拼起來得到的。理論上，如果假定一次投影的傅里葉變換形狀像一個(gè)切成薄片的“派”，我們可以簡單地把每個(gè)楔子插入適當(dāng)位置，以得到物體的一個(gè)二維傅里葉變換。不幸的是，每個(gè)投影傅里葉變換形狀類似在頻率空間的一個(gè)長條。如果簡單地計(jì)算所有投影傅里葉變換的和(假設(shè)在角度上等間隔)，中心區(qū)域被人為地增強(qiáng)，而外側(cè)區(qū)域數(shù)值不足。為了用條形區(qū)域估計(jì)“派”狀區(qū)域，我們可以給條形傅里葉變換乘以一個(gè)函數(shù)，該函數(shù)在靠近中心位置強(qiáng)度低，靠近邊緣時(shí)強(qiáng)度高。例如，可以將投影的傅里葉變換與該頻率處“派”狀楔子的寬度相乘。如果假設(shè)N個(gè)投影在180°內(nèi)均勻間隔，每個(gè)楔子的寬度在頻率是。權(quán)重函數(shù)的最終作用是加權(quán)長條的累加與“派”狀楔子的累加具有相同的“質(zhì)量”。[8]