基于光柵的X射線相位襯度成像圖像處理系統

魯展博 魏功祥

山東理工大學 山東 淄博 255049

一、研究背景及意義

傳統的X射線CT成像是基于樣品的吸收信息。如果樣品中含有大原子序數元素的樣品結構,如骨骼中的Ca元素,吸收信息明顯,利用傳統的基于吸收的CT成像可以很好的完成研究目標。但輕元素構成的組織結構器官對X射線具有相似的吸收截面,想要把病理結構從非病理結構中區分出來,用當前基于吸收襯度X射線成像的醫用X射線檢測系統就無能為力了。這些組織結構對X射線的吸收就非常弱了,無法利用傳統的CT得到足夠的襯度圖像。

研究發現,在硬X射線波段,輕元素引起的相位變化是其對硬X射線吸收的103-105倍,光學元件的相移截面要遠大于吸收截面。因此,要得到足夠襯度的相稱圖像就要比振幅成像方法所需的輻射劑量小得多。利用X射線的相位信息發展的相位襯度成像,特別適合對主要由輕元素組成的物體進行低輻射劑量觀測。

光波在物質中傳播,由兩個速度進行表征。一個是相速度,表示波節或波峰或任一相位狀態的傳播速度。群速度表示整個波包的傳播速度。光通過物質后,電子在光波的驅動下,圍繞原子核做受迫振蕩運動。

通常存在三種相位信號:(1)樣品像加速器或減速器,引起光波波陣面出現超前或滯后,與正相移或負相移信號對應;(2)像棱鏡物體,波陣面發生傾斜,導致光的折射,折射角與相位一階導數成正比;(3)樣品像透鏡,引起波陣面彎曲,導致光強的聚焦或發散,可以用相位的二階導數描述波陣面的彎曲程度。

由此發展了三種相位檢測方法:(1)利用干涉條紋探測相移;(2)利用角分辨元件探測折射角,獲得相位一階導數;(3)利用不同距離探測的光強獲得相位二階導數。

基于對軟組織成像的生物醫學研究要求,特別是腫瘤和癌細胞探測,需要一種能夠準確提取低Z元素組成的樣品結構信息的方法。由此發展了基于晶體干涉和基于傳播的同軸相襯等方法。

二、原理及依據

1、莫爾條紋的光強分布。周期為d,光柵線沿y方向的光柵透過率函數T(x,y)的Fourier級數表示:

an為第n階展開系數。光柵在距離為z的菲涅耳衍射區的衍射光場分布

在分數泰伯距離zT處:

對相位為π/2的相位光柵:在m為半整數,如1/2時,

此時,衍射光場分布表示為

2、物體相移引起的光場分布變化

其中,

為由于物體存在引起的微分相位.如果G1’和G2的光柵線的夾角為θ,則探測器所測得的莫爾條紋強度

圖1 光柵微分相襯成像原理圖

3、傅立葉變換法恢復微分相位

莫爾條紋光強分布可以簡寫為

其中,

對y坐一維傅里葉變換,得到其頻譜分布為

4 從微分相位通過希爾伯特濾波重構相位分布的理論

4.1 從光束偏折測量中進行斷層重建 從吸收或路徑長度測量等通常的線積分斷層重建算法是眾所周知的。斷層成像中最快捷和精確的算法是卷積后投影算法。該算法是重建的解析解,可以通過更精密的采樣間隔來更精確的接近完備解。這里討論如何擴展該算法到光束偏折測量中。

圖2 傅立葉切片定理原理圖

考慮物函數f(x,y)。通過物體z=z0平面,與x軸成θ角的普通投影可以表示為

其中x’-y’相對于x-y坐標系轉過θ。省略變量z,因為它對結果的推導(derivation)沒有影響。投影的一維傅里葉變換為

這里,用大寫字母表示頻域的函數和坐標,小寫字母則表示對應的空域坐標。可以看到,物函數f(x,y)可以從一組投影的傅里葉變換中重現出來

利用卷積理論,關于Y’的傅里葉變換可以表示為

該濾波函數如圖3。

圖3 卷積后投影濾波函數(a)和希爾伯特頻譜濾波函數(b)

卷積后投影算法可以擴展到光束偏折測量中。考慮折射率分布為

對小的傾角,可以近似的認為光束路徑是直的,分別考慮兩個正交的傾角。這樣,(x,y)面內的高斯光束的傾角可以近似為

對于圖4中的旋轉坐標系

圖4 傅立葉切片旋轉坐標系

由公式(16)可知

利用公式(26),(27)可得

其中,sgn(Y’)是符號函數

利用逆傅里葉變換,公式(28)可以寫作

在空間域中為

卷積后投影算法包括兩部分,卷積(或投影濾波)和后投影(或對所有投影進行積分)。對投影頻率成分的卷積步驟的影響是濾波函數的傅里葉變換。而對積分,后投影對連續連續無噪聲數據是沒有影響的。但是,對于離散投影,需要對后投影進行插值運算。對頻率成分的插值是各種各樣的,其影響可以通過插值濾波的選擇進行消除。由于噪聲的影響,后投影積分可以看作是一個平均過程,可以消除噪聲。頻率噪聲的平均效應依賴于噪聲的性質。不考慮特殊情況下的插值和后投影,對投影頻率分量的卷積后投影的純作用是由頻域中濾波函數K(Y’)給出的。比較常規和橫向梯度投影的濾波函數,可見對常規投影,濾波函數放大了高頻分量,降低了低頻分量,但對橫向梯度投影的濾波作用則是平均(flat)的。因為其他對斷層成像的解析形式與卷積后投影是等價的,這種濾波也適用于這些算法。這與橫向梯度測量(正交方向的微分)的投影(對一個方向的積分)和原始物體(對負頻成分變號)有相同的空間頻率權重是一致的。

這種對橫向梯度投影的不同濾波會對這種形式的測量方法產生重要的影響。一種是投影測量中包含噪聲的情況。如果投影中的樣品與物體中最高空間頻率的間隔足夠小,則在Nyquist頻率附近的信號就很小。所以這些信號的空間頻率成分會在低頻處放大,而在接近Nyquist頻率的部分接近于零。而噪聲源在這種情況下通常不是限帶的。如白噪聲和量子噪聲通常會平鋪到各頻率成分。對通常投影的卷積后投影重建,斜坡濾波減小了低頻成分的分布,SNR通常在這里是很高的,放大了高頻成分的分布,而在這些地方通常具有較小的SNR。為了減小高頻噪聲,可以選擇不同的窗函數,來降低對傳統的投影濾波函數的高頻響應。降低噪聲是以減小空間分辨率為代價的。對于橫向梯度投影的平場濾波,低頻成分的權重和高頻成分的權重是相同的。低頻成分對傳統的投影更有意義,高頻成分沒有增加。我們所有的從光束偏折投影的重建都沒有使用降低高頻分量的窗函數,因此噪聲一直是比較小的。對橫向梯度投影的頻域窗函數的簡單形式也是有用的。因為窗函數只對負頻分量產生π相移,對傅里葉重建技術活重建分析技術實現濾波相對也會簡單。

三、主要功能和實用情況

基于光柵的X射線相位襯度成像圖像處理系統包括四部分項目設計內容:

(1)光柵Talbot-Lau干涉儀原理展示;(2)動態光柵微分相位成像;(3)基于傅里葉變換的光柵相位CT成像;(4)掃描光柵相位CT成像。每部分都是一個獨立的GUI,可以通過按鈕控制項目內容的切換,并完成相應的項目內容。

1、光柵Talbot-Lau干涉儀原理展示

光柵微分相位襯度成像是利用光柵的自成像效應,又稱為泰伯(Talbot)效應。該系統設計目標是針對同步輻射光源條件下的實驗結果處理,所以在設計中主要針對上海光源X射線光柵相位成像系統。實驗儀器包括同步輻射光源、相位光柵和吸收光柵兩個光柵元件,樣品臺和X射線CCD組成。作為原理演示,其設計功能目標:

(1)通過光柵衍射的模擬,了解吸收光柵和相位光柵的衍射性質,特別是相位光柵的分數泰伯效應。由此了解光柵相襯成像的基本原理。

(2)計算不同光柵參數下的實驗條件。如根據不同的實際光柵參數,計算分數泰伯距離,為實驗提供參考。特別是在確定的光柵條件下,確定調整過程中需要用到的光柵距離等參數。

設計功能包括:

(1)根據參數,構畫光柵圖像。

A、光柵參數設置。包括設計光柵幅面大小參數設置,光柵周期參數,相位光柵和吸收光柵的選擇,相位光柵相位值的設置。能夠畫出吸收光柵和相位光柵。可以設置相位光柵的相位值。

B、構畫光柵二維構圖和三維構圖,選擇不同的著色模式和三維構圖的不同視角模式。

(2)觀察光柵頻譜結構。

(3)演示光柵在不同距離處的菲涅耳衍射圖樣。動態控制衍射距離,實時顯示衍射光場分布。可以實時構畫不同距離處的衍射光場的振幅、相位、實部和虛部等不同著色模式下的光場分布。演示光柵菲涅耳衍射的剖面結構;

(4)演示在上海光源利用光柵相襯成像實驗裝置所測量的部分莫爾條紋圖樣。

2、動態光柵微分相位成像實驗數據處理系統

該部分的設計功能主要進行動態樣品的微分相位成像數據處理。該方法是基于光柵莫爾條紋,利用傅里葉變換算法實現樣品微分相位信息的提取。直接對樣品的微分相位進行成像。

本算法的難點在于利用光柵莫爾條紋對樣品相位信息的調制,怎樣通過適當的傅里葉變換算法,快速提取光柵調制的微分相位信息。要求能夠處理CCD的暗電流噪聲,光源的背景噪聲,光柵噪聲等系統噪聲,并能夠準確地提取樣品的微分相位頻譜,最大限度地復原樣品的微分相位信息。

項目組成目標:讀入光柵微分相位成像數據,自動進行系統噪聲處理。

功能目標包括:

(1)輸入基本實驗參數,并自動形成保存文檔;

(2)選擇實驗數據保存路徑,自動讀入實驗數據;

(3)實驗數據自動處理,并實時顯示;(由于部分實驗數據處理方法還處于保密階段,數據處理算法不能詳細公開)。

(4)自動保存實驗數據,制作動畫演示文稿。

3、基于傅里葉變換的光柵相位CT成像

本部分的設計功能是利用傅里葉變換算法準確地提取樣品的微分相位信息,并利用該方法對樣品進行CT重現。

傳統的CT重構是基于傅里葉切片定理,利用采集的不同角度的樣品的吸收(衰減)圖像數據,進行斷層切片重構。傳統的處理方法是利用濾波反投影處理。但對于通過處理所得到的光柵微分相位信息進行相位信息重構,需要利用希爾伯特變換濾波反投影算法,從微分相位信息直接重構相位信息。其關鍵技術在于希爾伯特變換濾波。

利用光柵微分相位數據,其主要優勢在于可以從一組數據中直接進行提取衰減信息,相位信息和散射信息。

本部分的項目目標是:

(1)自動進行數據的背景噪聲處理,得到投影圖;

(2)讀取投影文件總數,以及180°投影的平場數據采集周期,自動計算實驗數據,得到重建的微分相位數據組,振幅數據組和吸收信息數據組。并保存實驗復振幅數據,便于后續過程計算。

(3)sino構圖。根據傅立葉切片定理,由不同角度的投影數據,組合不同切片的投影數據組。利用Radon變換,得到切片數據組slice。

(4)CT重構。

A:由重建的微分相位成像的切片數據組進行希爾伯特變換,得到相位襯度CT數據組。

B:根據重建的吸收像的切片數據組進行濾波反投影,得到吸收CT數據組。

C:根據重建的散射像的切片數據組進行濾波反投影,得到散射CT數據組。

(5)根據CT數據組,進行三維構圖。

4、掃描光柵相位CT成像

目前,最精確的光柵相位襯度成像算法是在一個光柵周期中得到幾幅樣品的莫爾條紋,利用各像素的振蕩曲線,準確恢復該像素點的微分相位信息。其主要缺點是需要在一個光柵周期中進行掃描,所以CT掃描速度無法提高,CT設備也更為復雜。其優點是計算精度較高。

本部分的設計功能是循環讀取光柵掃描振蕩圖樣,從每個周期的掃描圖樣中提取各像素的振蕩曲線,并解出該像素的微分相位數據;將解出的微分相位數據作為投影數據進行CT重構。

4、基于光柵的X射線相位襯度成像圖像處理系統展示

1、成像及圖像處理系統

在上海同步輻射光源搭建得光柵成像系統如圖5.1所示,系統主要由樣品旋轉臺、光柵位移臺、相位光柵、吸收光柵、X射線探測器以及PC控制器這幾個裝置構成。其中相位光柵的規格為π/2的相位調制光柵,吸收光柵則使用的規格為柵格周期2.4μm,鍍金厚度為60μm的光柵,兩光柵沿橫向進行掃描。

圖5 上海光源光柵相稱成像系統

圖6展示的便是動態光柵微分相位成像實驗數據處理系統首頁界面,功能主要進行動態樣品的微分相位成像數據處理。該方法是基于光柵莫爾條紋以及利用傅里葉變換算法實現樣品微分相位信息的提取。直接對樣品的微分相位進行成像。利用傅里葉變換算法準確地提取樣品的微分相位信息,并利用該方法對樣品進行CT重現。

圖6 基于光柵的X射線微分相位襯度成像(英/中文)首頁

2、光柵的Talbot效應演示

功能說明:

主要包括4個控制模塊:

(1)參數設置。通過設置M,N兩個參數,調整光柵尺度。Period調整光柵周期像素數,光柵都采用羅切型光柵,占空比1:1。Tilt angle調整光柵線傾斜角度。通過選擇phase和absorb選擇光柵類型是相位型,還是吸收型光柵。如果選擇相位型光柵,可以設置光柵的相位值參數,通過下拉列表框,給出1、π/2,π和random四種相位選擇值。光柵參數設置完后,可以通過按鈕”2D draw”和’3D draw”進行二維和三維構圖,并可以通過colormap選擇調色板。如果選擇復選框rotate 3D,可以對三維圖像進行三維旋轉。

(2)頻譜分析。光柵頻譜結構可以通過按鈕Frequency進行計算,本算法是對光柵進行FFT進行計算的,顯示采用二維振幅顯示。由于頻譜數值范圍較大,在畫圖過程中采用對數歸一化方式顯示。為了定量的觀察頻譜結構,可以通過按鈕”Profile”,選擇觀察頻譜中心部分的剖面結構,通過scale調整顯示的振幅范圍,并通過”position shift”滑塊調整顯示坐標位置。如果中心頻譜影響頻譜細微結構的觀察,可以通過”delete zero frequency”切掉零級頻譜進行觀察。

圖7 光柵泰伯效應演示操作頁面

(3)光柵衍射。

光柵衍射部分的主要功能:

1)根據上海光源的能量計算波長,計算公式為

λ=1.24/E

波長單位nm,Energy單位keV。

2)計算光柵模擬實驗中的像素尺寸。光柵周期在功能(1)中設置,像素實際大小dx=實際光柵周期d/模擬光柵周期像素數pn。模擬光柵周期像素數為功能(1)中的period參數,實際光柵周期為“actual grating period(um)”參數,單位um。像素實際大小為”pixel scale”,單位為um。光柵一級泰伯效應距離計算公式為

Td=2d2/λ

單位一般采用mm表示。

3)在光柵相襯實驗中,通常利用相位光柵的分數泰伯距離產生振幅型光柵自成像,需要觀察分數泰伯距離處的光柵自成像效果。所以在本方案中增加了分數泰伯距離的計算。可以通過點選fraction選擇分數值,并計算得到分數泰伯距離。

4)光柵衍射。通過拖動滑塊調整不同衍射距離,本算法采用自適應標量衍射理論編寫衍射函數[]。

默認顯示衍射圖樣的振幅,可以通過選擇“show item”,選擇顯示衍射復振幅數據的振幅、相位、實部或者虛部。

圖8 光柵衍射參數設置與計算

圖9中所示分數值為0時的光柵振幅像,也就是光柵本身。分數值為1的Talbot自成像。可以看到中間區域能夠很好的自成像,但在邊緣處會出現較大偏差,這是由于光柵邊緣部分不具備周期性所致。分數值為1/2處的自成像。由于光柵占空比為1:1,在1/2Talbot距離處,也可以得到較好的光柵自成像。光柵在2倍Talbot距離處的自成像。自成像的偏差主要是由于模擬中光柵的有限大小所造成的。

圖9 (a)光柵分布的振幅像(b)分數值為1的吸收型Talbot自成像(c)吸收型光柵在1/2Talbot距離處的自成像又稱分數泰伯效應(d)吸收型光柵在2倍Talbot距離處的自成像

5)不同載頻的莫爾條紋和實驗記錄結果

光柵在2倍泰伯距離處的自成像。自成像的偏差主要是由于模擬中光柵的有限大小所造成的。

(4)衍射剖面

光柵剖面部分的功能是給出光柵衍射過程中,光柵某個位置的衍射演變過程。

圖10 不同載頻的莫爾條紋和實驗記錄結果

3、動態光柵微分相位襯度成像處理步驟

該功能主要是進行動態樣品的微分相位成像數據處理。方法是基于光柵莫爾條紋。如圖12所示的便是動態光柵微分相位襯度成像處理過程整個的一個流程,主要包括暗電流噪聲的矯正、背景噪聲提取、物體信息提取、物體信息再現的環節。

圖11 相位為π和π/2的相位光柵在不同距離處的光柵衍射剖面

圖12 動態光柵微分相位襯度成像處理流程

圖13 動態光柵微分相位襯度成像的處理界面

其系統操作界面如圖12所示,這里以活體螞蟻為樣本進行了動態光柵微分相位襯度成像其圖像處理結果如圖15,圖14所示的分別為背景矯正的過程和頻譜濾波截取的過程。

圖14 (a)背景矯正界面示意圖(b)頻譜濾波截面演示示意圖

圖15 (a)再現相位像(b)再現振幅像(c)再現吸收像

如圖16所示,可以通過按鈕控制圖像處理進程,實現停止,暫停,重做,退出等。并自動保存處理圖像,可根據需要最后得到活體螞蟻動態微分相位信息視頻。4、基于傅里葉變換的光柵相位CT成像

圖16 自動批處理過程

(1)界面

圖17 光柵微分相位襯度CT成像操作界面

(2)操作與處理結果

圖18 再現結果。(a)經過矯正的莫爾條紋,(b)再現吸收像,(c)再現微分相位像,(d)散射像

五、創新點和特色

1、首次實現了利用傅立葉變換方法,從一幅記錄的莫爾條紋圖像中恢復樣品的微分相位信息。因為不需要掃描,可以實現相位的快速CT成像,有利于減小生物樣品的輻照劑量。由于實驗過程中,光源和光柵的非理想化,都會對實驗結果引入很強的噪聲信息,對相位信息的提取造成非常大的干擾。噪聲來源包括暗電流噪聲,平場噪聲,光柵噪聲。本項目的主要難點是分析這些噪聲來源,并通過適當的處理步驟消除噪聲干擾。

2、能夠方便的選擇合適的條紋載頻和頻譜濾波窗口。頻譜濾波是一個非常敏感的操作,而莫爾條紋的頻譜結構要比一般的全息圖的頻譜復雜很多,如何從莫爾條紋的頻譜信息中獲取所需要的頻譜結構,是該方法的關鍵。濾波窗口太大,會將一些不需要的頻譜結構保留下來,重疊在恢復的相位信息中,造成相位信息的不準確。濾波窗口太小,會損失樣品的高頻信息,使再現的樣品信息不完整,甚至會出現較嚴重的震蕩效應。

3、首次利用希爾伯特濾波反投影算法,從恢復的微分相位信息獲取斷層相位信息。與傳統的濾波反投影算法不同,濾波反投影適用于吸收信息的再現。而我們所得到是微分相位信息,如果對每一幅微分相位圖像進行相位恢復,CT重建過程會變得非常慢。因此,我們提出在濾波過程中,直接采用希爾伯特反投影算法,可以直接得到相位斷層信息,處理速度明顯加快。

注 釋

[1]Gregory W.Faris and Robert L.Byer,Three-dimensional beam-deflection optical tomography of a supersonic jet.APPL.OPT.(1988)27(24),5202

[2]國承山等,光衍射數值模擬中不同抽樣方法的適用性分析。光學學報,2008,28(3)442-446”