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

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

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

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

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