通常,，造影劑注入前的DSA血管圖像被稱為蒙片(Mask Image)，而同一部位在注入造影劑后的X射線圖像被稱為活片(Live Image)。在假設(shè)血管周圍組織,、成像條件完全不變的情況下,，將蒙片與活片相減應(yīng)得到僅包含血管的圖像。但在臨床應(yīng)用中,，這種理論假設(shè)總是不成立的,。在成像的同時(shí)，病人會(huì)有各種運(yùn)動(dòng),，其中有些運(yùn)動(dòng)是不可避免的,，如呼吸、肌肉運(yùn)動(dòng),、心臟運(yùn)動(dòng)以及在造影劑注入時(shí)由于病人的自然反應(yīng)而引起的局部運(yùn)動(dòng)等,。這些運(yùn)動(dòng)會(huì)使活片與蒙片之間存在并非由造影劑造成的差異，從而使減影后的圖像出現(xiàn)偽影,，妨礙醫(yī)生的正確診斷,。本文提出了一種自動(dòng)消除數(shù)字減影血管造影DSA(Digital Subtraction Angiography)減影圖像中偽影的方法。通過(guò)實(shí)驗(yàn)證實(shí),，本方法比其他類似的方法^[1]更快,、更有效。
1 方法原理
　　DSA序列圖像的配準(zhǔn)目的是通過(guò)圖像處理和圖像分析手段找到一種幾何變換方法能消除序列圖像中任意二幅之間因?yàn)椴∪说倪\(yùn)動(dòng)帶來(lái)的偽影,。本文選擇對(duì)蒙片作相應(yīng)的變形(warping),，使其在圖像相減之前與造影圖像相匹配,。
　　DSA圖像的實(shí)質(zhì)就是X射線投影圖像,。根據(jù)不同物質(zhì)對(duì)X射線的吸收系數(shù)的不同，在三維空間把沿著射線方向的不同物體的吸收系數(shù)疊加,，形成二維投影圖像,。因此，DSA圖像中的運(yùn)動(dòng)偽影是由三維空間運(yùn)動(dòng)造成的,。
2 實(shí)現(xiàn)過(guò)程
　　配準(zhǔn)的具體過(guò)程為：首先,，對(duì)于一幅圖像中的某一點(diǎn)，根據(jù)一定的相似度函數(shù)找到另一幅圖像中的對(duì)應(yīng)點(diǎn),，并求出" title="求出">求出相對(duì)位移,；然后對(duì)其中一幅圖像(蒙片)進(jìn)行變形，以達(dá)到與另一幅圖像(活片)配準(zhǔn)的目的,。但如果分別計(jì)算每一點(diǎn)的位移,，則計(jì)算量太大，不能滿足臨床要求,。因此可以先求少數(shù)特征點(diǎn)(控制點(diǎn))的位移,，再通過(guò)插值方法得到所有點(diǎn)的位移，以此得到快速配準(zhǔn)的方法。下面介紹具體算法,。
2.1 找控制點(diǎn)
　　控制點(diǎn)的選取有多種方法,，如用人工方法選取感興趣區(qū)域中的某些點(diǎn)，或者在圖像中用規(guī)則的網(wǎng)格選取等,。經(jīng)觀察,，大多數(shù)運(yùn)動(dòng)偽影出現(xiàn)在圖像的邊緣點(diǎn)上，因此可以認(rèn)為邊緣點(diǎn)的匹配比均勻區(qū)域的匹配更可靠,，選擇用邊緣點(diǎn)作為控制點(diǎn),。
　　通過(guò)邊緣檢測(cè)" title="邊緣檢測(cè)">邊緣檢測(cè)算子可以找到圖像的邊緣點(diǎn)。邊緣檢測(cè)方法很多,，本文選用Canny邊緣檢測(cè),。它是一種比較新的檢測(cè)算子^[2]，具有很好的邊緣檢測(cè)性能,，正得到越來(lái)越廣泛的應(yīng)用,。Canny邊緣檢測(cè)法利用高斯函數(shù)的一階微分，能在噪聲抑制和邊緣檢測(cè)之間取得較好的平衡,。由于沒(méi)必要把所有用Canny算子檢測(cè)出的點(diǎn)都作為控制點(diǎn),，所以在Canny邊緣點(diǎn)計(jì)算出該點(diǎn)的梯度值后，根據(jù)閾值來(lái)判斷該點(diǎn)是否作為控制點(diǎn),。
　　為使控制點(diǎn)的分布更合理,，利用比例因子對(duì)控制點(diǎn)間的距離作規(guī)定：每2個(gè)控制點(diǎn)間的最小距離為D_min=λ_minM。同時(shí)為避免在很大區(qū)域內(nèi)沒(méi)有控制點(diǎn),，也規(guī)定控制點(diǎn)間的最大" title="最大">最大距離：D_max=λ_maxM,。通常在DSA圖像中，有用的信息都顯示在一個(gè)圓形區(qū)域內(nèi),，區(qū)域之外的像素值為零,。這樣的圓形區(qū)域被稱為“曝光區(qū)域”，其定義如下：
　　R=(x-x_c)²+(y-y_c)²≤r²,，X_min≤x≤X_max,，Y_min≤y≤Y_max
　　其中，(x_c,，y_c)為影像增強(qiáng)器的中心點(diǎn),，通過(guò)適當(dāng)校正可與圖像的中心點(diǎn)重合。r為圓形區(qū)域的半徑,，而X_min,、X_max、Y_min,、Y_max分別是圖像的左,、右和上,、下邊界。既然所有有用信息都在圓形區(qū)域內(nèi),，因此控制點(diǎn)也應(yīng)該在R之內(nèi)尋找,。
　　圖1為DSA蒙片及其控制點(diǎn)圖像，其中圖1(a)為一腦部DSA正面圖像的蒙片,，圖1(b)是從圖1(a)中得到的控制點(diǎn)圖像,。由圖可見(jiàn)，控制點(diǎn)都位于邊緣上,，且保持著一定的均勻度,。

2.2 選擇相似性測(cè)度
　　本文采用模版匹配法求蒙片與活片中像素的相對(duì)位移：=(d_x，d_y),。模版匹配法利用相似性測(cè)度來(lái)求一幅圖像I₀(x,，y，t₀)中的某一點(diǎn)P(x,，y)相對(duì)于另一幅圖像I₁(x,，y，t₁)(t₁＞t₀)的位移,，以P點(diǎn)為中心,，以I₀中P點(diǎn)周圍W×W鄰域內(nèi)的點(diǎn)為一窗口，在I₁中求出這一鄰域移動(dòng)后的窗口位置,，并利用適當(dāng)?shù)南嗨菩詼y(cè)度,，求出使相似性測(cè)度達(dá)到最大的位移，即為要求的最佳位移,。
　　本文選用基于圖像直方圖的相似性測(cè)度^[3],，當(dāng)兩幅圖達(dá)到最佳匹配時(shí)，兩幅圖像之差的直方圖會(huì)集中在很窄的一段灰度范圍內(nèi),，且形成明顯的峰值,，而當(dāng)兩幅圖像不是很匹配時(shí),，直方圖的集中程度會(huì)降低,；更重要的是，直方圖的這種性質(zhì)不會(huì)受血管的影響,，DSA造影圖像及其在不同配準(zhǔn)情況下的直方圖比較如圖2所示,。其中圖2(a)是DSA造影圖像，其中的白框表示被選取的區(qū)域,，左上方的白框不包括血管,，而中間的白框包括血管；圖2(b)和圖2(c)分別表示不包括血管時(shí),，最佳和非最佳配準(zhǔn)下差異圖像的直方圖,；圖2(d)和圖2(e)分別表示包括血管時(shí)，最佳和非最佳配準(zhǔn)下差異圖像的直方圖。

　　假設(shè)待配準(zhǔn)的兩幅圖像分別為I₀(x,，y,，t₀)和I₁(x，y,，t₁)(t₁＞t₀),，求I₀中某一點(diǎn)P₀的位移。在P₀周圍選取一個(gè)小窗口W₀,，該窗以P₀為中心,，大小為W×W個(gè)像素，對(duì)于每一個(gè)位移在I₁中找到相應(yīng)的相同大小的窗口W1,，將這兩個(gè)窗口相減,，得到差異圖像Wd(W×W)，針對(duì)Wd計(jì)算相似性測(cè)度,，有研究人員提出基于直方圖的判定標(biāo)準(zhǔn)^[3]：
　　

　　其中[δ_min,，δ_max]?奐Z是差異圖像中可能的灰度范圍，h(δ)是差異圖像的歸一化直方圖,，并假設(shè)f是嚴(yán)格凸函數(shù)或嚴(yán)格凹函數(shù),，但可微。相關(guān)文獻(xiàn)提出了多種基于直方圖的f函數(shù)形式^[4],，并指出當(dāng)f(x)=x²時(shí)計(jì)算簡(jiǎn)單且精度高,。這樣，式(1)應(yīng)寫成：
　　

2.3 亞像素精度
　　某些研究者指出,，即使是亞像素級(jí)的移動(dòng)也會(huì)造成兩幅圖像間的差異很大^[5],，因此在計(jì)算像素位移時(shí)應(yīng)設(shè)法達(dá)到亞像素精度，而這要用到圖像插值技術(shù),。本方法通過(guò)先對(duì)圖像插值和重采樣,，然后對(duì)插值后的圖像求相似性測(cè)度來(lái)達(dá)到亞像素精度。期望達(dá)到的亞像素精度為0.1 pixel,。
　　圖像插值在整個(gè)配準(zhǔn)算法中很重要,。它的精度直接影響著位移的精度；它的計(jì)算速度影響著整個(gè)算法的計(jì)算時(shí)間,。如何選擇一種既精確又快速的插值方法成為此配準(zhǔn)算法成敗的關(guān)鍵,。本文采用三次多項(xiàng)式插值法對(duì)原圖像I₀進(jìn)行插值。為得到理想的效果,，本文采用4×4的卷積核,。
2.4 最優(yōu)化搜索
　　在確定了相似性測(cè)度為直方圖能量E()后，要求最佳位移使相似性測(cè)度達(dá)到最大,。本問(wèn)題屬于無(wú)約束優(yōu)化問(wèn)題,，針對(duì)該問(wèn)題,，本文采用直接搜索法。最常用的是鮑威爾(Powell)法^[6]和單純形法,。這里采用鮑威爾法,，也稱為方向加速法。
　　搜索的速度受計(jì)算位移時(shí)選用的窗W×W大小影響顯著,，這是個(gè)矛盾的統(tǒng)一體：一方面為了減少搜索時(shí)間,，希望采用較小的窗；另一方面,，較小的窗將使相似性測(cè)度產(chǎn)生許多局部極大點(diǎn),，導(dǎo)致搜索結(jié)果不可靠。因此在選擇窗的大小時(shí)應(yīng)該慎重考慮,，本文選用的窗口大小為51×51個(gè)像素,。
2.5 位移插值
　　求出所有控制點(diǎn)的最佳位移后，進(jìn)而需要求出圖像I₀中所有點(diǎn)相對(duì)I₁的位移,，這要利用線性插值來(lái)估算出所有點(diǎn)的位移,。本文采用Delaunay三角插值法^[7]，它保證在每一個(gè)三角形內(nèi)的最小角度盡可能地最大化,，且保證了這種三角形分法的惟一性,。
　　獲得圖像I₀中所有點(diǎn)的位移后，即可對(duì)圖像I₀進(jìn)行校正,。因?yàn)槲灰?IMG src="http://files.chinaaet.com/old/uploadfiles/jishu/jslw/20080701034728234.gif" border=0>的單位是0.1像素,，則圖像上一點(diǎn)根據(jù)移動(dòng)后極可能落入整數(shù)值的網(wǎng)格之間，這里需要用插值法(雙線性插值)得到非整數(shù)像素位置上的像素值,。
3 實(shí)驗(yàn)結(jié)果及討論
　　拍攝后前位的DSA時(shí)間序列圖像共39幅,，每幅圖像的大小為512×512個(gè)像素，選取其中的第1幅(蒙片)和第28幅(活片)作DSA圖像的配準(zhǔn)研究,，其中的蒙片如圖1(a),。
　　控制點(diǎn)提取之后(見(jiàn)圖1(b))，根據(jù)模板匹配法計(jì)算每個(gè)控制點(diǎn)的最佳位移,，亦即使差異圖像的直方圖能量達(dá)到最大時(shí)的位移,。此時(shí)，對(duì)窗口大小W的選擇很重要,。若采用鮑威爾最優(yōu)化方法搜尋最佳位移,，則要求E(d)足夠光滑以至于只有一個(gè)極值,。窗口太小將增加E(d)的局部極小值,，從而影響到最優(yōu)解的精度。經(jīng)多次嘗試和比較,，選用的窗口大小為51×51個(gè)像素,。