本發(fā)明涉及一種醫(yī)學圖像處理，具體涉及圖像的增強或復原，特別是涉及一種圓軌道下不完全ct投影數(shù)據(jù)的恢復方法。
背景技術(shù)：
：錐形束計算機斷層掃描是一種可以快速且低輻射劑量獲得人體的三維解剖圖像的醫(yī)學成像方法。數(shù)字平面探測器的出現(xiàn)使得真正的錐形束計算機斷層掃描(cone-beamcomputedtomography,cbct)在臨床應用成為可能。由于cbct采用錐形束掃描結(jié)構(gòu)，所采集的投影數(shù)據(jù)具有高度冗余性，想要利用數(shù)據(jù)的冗余性來達到高效率的目的就必須滿足一致性條件(consistencyconditions)。約翰于1938年從數(shù)學的角度推導出了一個超雙曲線偏微分方程(ultrahyperbolicpartialdifferentialequation,pde)作為線積分的一致性條件，故此方程稱為約翰方程。在70年代，隨著ct的出現(xiàn)和興起，約翰方程(john'sequation)被廣泛地應用在ct領域并進行了更深入地推導和延伸。時至今日，已經(jīng)有許多學者提出了多種基于約翰方程的一致性條件，如patch提出一種基于第三代螺旋ct的一致性條件(sarahk.patch,computationofunmeasuredthird-generationvctviewsfrommeasuredviews,2002,ieeetrans.med.imaging)，在此文獻中，作者通過使用第三代螺旋ct參數(shù)來替換基于笛卡爾坐標系的約翰方程中的參數(shù)，并在此基礎上進行了一系列地深入推導，將二階偏微分方程轉(zhuǎn)換成了一組一階常微分方程(first-orderordinarydifferentialequation)，對此方程進行傅里葉變換，然后應用其頻率域的形式可以從獲得的投影中計算出未測量的投影數(shù)據(jù)。根據(jù)此理論，文獻中使用“旋轉(zhuǎn)+螺旋+旋轉(zhuǎn)”的射線源軌道，螺距為0.128m，每轉(zhuǎn)一圈采集984個投影，然后針對獲得的投影數(shù)據(jù)使用推導方程對應的傅立葉變換形式進行外插，從而獲得未測量的投影數(shù)據(jù)。但是patch推導出的一致性條件是基于螺旋ct的幾何結(jié)構(gòu)推導而成，螺旋結(jié)構(gòu)要求變量z沿著縱軸在不停地變化，這會在公式中引入一個變化量，使得公式計算復雜，解決問題困難。文獻中由于采用“旋轉(zhuǎn)+螺旋+旋轉(zhuǎn)”掃描方式，需要在“螺旋”這個掃描階段前后分別進行工作臺移動速度的增加和減小，以保證在這三個掃描階段數(shù)據(jù)采集的快速轉(zhuǎn)換和平滑過渡，但是在實際應用中，一般螺旋ct的工作臺移動速度都是恒定的，因此，文獻中的一致性條件很難具體應用。此外，haoyan等人在使用移動的波束停止陣列(beamstoparray,bsa)來進行散射校正時采用了patch基于三代螺旋ct的一致性條件(haoyan,xuanqinmou,shaojietang,qiongxuandmariazankl,projectionmovingbeamstoparray,2010,phys.med.biol.)。波束停止陣列是一個遮擋板，上面有規(guī)律地排列著遮擋點。x射線經(jīng)過遮擋板照射到物體上，其中經(jīng)過遮擋板上遮擋點的x射線會被屏蔽掉。如果在探測器遮擋點對應位置出現(xiàn)數(shù)值，則這個數(shù)值是散射值。由于使用了bsa，得到的投影數(shù)據(jù)是不完整的，因此需要使用一致性條件對其進行恢復。此文獻采用patch基于三代螺旋ct的推導公式[patch(2002)]，并結(jié)合空間插值方法，提出一種pc_si方法用于散射校正領域中。但如前面提到的那樣，由于螺旋ct的z軸在不停地變化，導致此公式很難具體應用，所以此文獻參照defrise等人的做法，用探測器的縱坐標來近似代替z坐標，但是這種近似做法，會導致精度的下降。從幾何結(jié)構(gòu)上來看，只需要令z＝0，螺旋ct就可以退化成圓軌道cbct，但是文獻中的公式是對z求導，當z為常數(shù)時導數(shù)沒有意義，所以此文獻中所用的方法并不適用于圓軌道cbct。因此，針對現(xiàn)有技術(shù)不足，提供一種圓軌道下不完全ct投影數(shù)據(jù)的恢復方法以克服現(xiàn)有技術(shù)不足甚為必要。技術(shù)實現(xiàn)要素：本發(fā)明所要解決的技術(shù)問題是提供一種圓軌道下不完全ct投影數(shù)據(jù)的恢復方法，和傳統(tǒng)的空間插值方法相比，本發(fā)明的方法能夠很好地恢復丟失信息中的高頻成分，進而大幅度減少圖像中的條形偽影，顯著提高圖像質(zhì)量。在圓軌道下使用patch的方法需要先使用探測器的縱坐標來近似代替z坐標，而本方法采用的一致性條件是基于圓軌道cbct幾何結(jié)構(gòu)推導得到的，能夠避免因為近似z坐標而導致的精度下降，并且patch的方法中z軸在不停變化，使得此方法具體應用困難。與之相反，本發(fā)明方法只是一個簡單的一階求導公式，實現(xiàn)簡單，具有可應用性。本發(fā)明采用如下技術(shù)方案：一種圓軌道下不完全ct投影數(shù)據(jù)的恢復方法，通過如下步驟進行：第一步，獲得壞像素點位置信息；第二步，對每個投影中的壞像素點使用水平一維三次樣條插值方法進行線性插值，以恢復壞像素點的低頻信息，經(jīng)過此步驟得到每個角度下的初始化投影；第三步，令n＝1，進入第四步，第四步，令m＝n，m為當前迭代次數(shù)；將每個角度下的初始化投影進行二維傅里葉變換，獲得每個角度下的初始化投影在頻率域中的頻域數(shù)據(jù)；第五步，分別將前后兩個相鄰投影的頻率域形式帶入到一致性條件對應的公式中計算投影結(jié)果；第六步，對加權(quán)后的投影進行傅里葉反變換，得到第m次迭代終值；第七步，判斷當前迭代次數(shù)m是否大于等于迭代終止次數(shù)s，如果是，則進入第九步；否則進入第八步；第八步，將第m次迭代終值作為初始化投影，令n＝n+1，進入第三步；第九步；以第m次迭代終值為最終恢復的數(shù)據(jù)結(jié)果。優(yōu)選的，上述的圓軌道下不完全ct投影數(shù)據(jù)的恢復方法，第五步中使用的一致性條件，具體形式如下：其中，其中ρ是射線源到旋轉(zhuǎn)中心的距離，d是旋轉(zhuǎn)中心到探測器的距離，(α1,α2)是平板探測器坐標，θ是方位角。優(yōu)選的，上述第五步具體是使用相同的權(quán)重因子ω對前后兩個投影結(jié)果進行加權(quán)，得到加權(quán)更新后的投影結(jié)果。優(yōu)選的，上述的圓軌道下不完全ct投影數(shù)據(jù)的恢復方法，第五步，具體使用一致性條件即公式(28)對每一個初始化投影的頻域數(shù)據(jù)進行恢復，公式(28)的具體形式如下：u*(θ+dθ)＝u*(θ)+dθ·u,k2≠0……公式(28)；其中，u*(θ+dθ)和u*(θ)分別表示θ+dθ和θ角度下初始化投影的頻域數(shù)據(jù)，u*(θ)為已知量，u*(θ+dθ)是要求的未知量，dθ代表相鄰兩個角度的間隔，k1、k2是平板探測器的頻率域坐標，i代表復數(shù)的虛部；公式(28)的具體使用方法如下：首先將u*(θ+dθ)代入到公式(28)中計算u*(θ)：u*(θ)＝u*(θ+dθ)-dθ·u；然后用u*(θ-dθ)代替公式(28)中的u*(θ)，u*(θ)代替公式(28)中的u*(θ+dθ)計算得到：u*(θ)＝u*(θ-dθ)+dθ·u；由于選取前后兩個相鄰投影頻域數(shù)據(jù)的任一個都可以實現(xiàn)投影數(shù)據(jù)恢復，為了平衡恢復效果，使用權(quán)重因子ω，ω∈[0,1]，對上述兩個公式的計算結(jié)果進行加權(quán)平均，得到加權(quán)后的結(jié)果ur*(θ)，即:ur*(θ)＝ω(u*(θ-dθ)+dθ·u)+(1-ω)(u*(θ+dθ)-dθ·u),k2≠0；由于公式(28)只有在k2≠0時才有效，所以對頻域數(shù)據(jù)u*(θ)中k2＝0的部分，直接使用θ角度下初始化投影頻域數(shù)據(jù)中的相應部分來補償；第六步，具體是對ur*(θ)進行傅里葉反變換，得到第m次迭代終值ur(θ)；ur(θ)＝f-1(ur*(θ))……(1)；其中公式(1)中ur(θ)表示θ角度下加權(quán)后的頻域數(shù)據(jù)ur*(θ)的反傅里葉變換結(jié)果，f-1代表傅里葉反變換；第七步，具體是判斷當前迭代次數(shù)m是否大于等于迭代終止次數(shù)s，如果是，則進入第九步；否則進入第八步；第八步，具體是將第m次迭代終值ur(θ)作為初始化投影，令n＝n+1，進入第三步；第九步，具體是以第m次迭代終值ur(θ)為最終恢復的數(shù)據(jù)結(jié)果。優(yōu)選的，上述的圓軌道下不完全ct投影數(shù)據(jù)的恢復方法，第一步，獲得壞像素點位置信息具體是將實際采集的投影數(shù)據(jù)轉(zhuǎn)換成對應的弦圖，即：弦圖的x軸為投影的像素點位置，y軸為該像素位置對應的投影角度位置；對某個探測器單元存在損壞的條件下，則在弦圖中對應的像素位置上會出現(xiàn)一條豎線，根據(jù)此理論在弦圖中獲得壞像素點的位置信息，并記錄壞像素點位置。其中本方法中應用的公式(28)推導過程如下：由于本發(fā)明提出的一致性條件是基于約翰方程，所以先對約翰方程進行簡單介紹和分析。fjohn于1938年在文獻《theultrahyperbolicdifferentialequationwithfourindependentvariables》中提出一個超雙曲線二階偏微分方程，作為線積分的一致性條件，具體形式如下：公式中i,j為下標；u是函數(shù)f穿過ε和η兩點的線積分：u(ε；η)＝∫f(ε+t(η-ε))dt(3)；以下為圓軌道下基于約翰方程的一致性條件推導步驟：首先采用圓軌道cbct坐標參數(shù)對ε和η進行參數(shù)化：其中ρ是射線源到旋轉(zhuǎn)中心的距離，d是旋轉(zhuǎn)中心到探測器的距離，(α1,α2)是平板探測器坐標，θ是方位角。分別求ε和η的每個分量求偏導：為書寫方便，定義如下算子：把公式(6)-(11)代入到公式(2)中得到：由公式(13)-(16)對h1，h2,的定義，代入到公式(17)-(19)中，得到：從上述公式中，得出：對上面公式(21)移項合并同類項：求解公式(23)，最終推導結(jié)果如下：公式(24)即為我們推導的圓軌道下基于約翰方程的一致性條件。對公式(24)進行傅里葉變換，得到對應的頻率域的表達形式如下：在上述公式中，uθ*表示u對θ求導后的傅里葉變換，(k1,k2)是(α1,α2)的對應頻率域形式。為了方便表達和應用，將公式(25)右邊式子用u代替：考慮到求導公式的一般定義，可以將公式(25)改下成如下形式：即：u*(θ+dθ)＝u*(θ)+dθ·u,k2≠0(28)；公式(28)是本發(fā)明改進后的一致性條件，本發(fā)明使用公式(28)進行投影數(shù)據(jù)的恢復。傳統(tǒng)的空間插值技術(shù)使用同一個角度的投影的其他完整像素點數(shù)據(jù)來恢復丟失的數(shù)據(jù)，所使用重建信息只局限于丟失數(shù)據(jù)所在角度的投影，所以插值方法只能恢復丟失數(shù)據(jù)中的低頻信息，對丟失的高頻信息基本沒有恢復。本發(fā)明方法應用的一致性條件，即不僅基于同一個投影中的數(shù)據(jù)來恢復丟失數(shù)據(jù)，還利用相鄰角度投影中的對應位置數(shù)據(jù)信息。因此本發(fā)采用一致性條件可以很好地恢復丟失數(shù)據(jù)中的高頻信息，能基本消除使用插值技術(shù)獲得的重建圖像中的條形偽影，顯著提高圖像質(zhì)量。附圖說明圖1為實例中圓軌道cbct掃描幾何示意圖及參數(shù)示意圖，圖中，ρ是射線源到旋轉(zhuǎn)中心的距離，d是旋轉(zhuǎn)中心到探測器的距離，(α1,α2)是平板探測器坐標，θ是方位角。圖2是投影數(shù)據(jù)對應的弦圖，橫坐標x軸為像素點位置，縱坐標y軸為投影角度。圖3為本發(fā)明的數(shù)據(jù)恢復方法流程圖。圖4為三種方法的仿真結(jié)果示意圖，圖4(a)為理想圖像，圖4(b)為未經(jīng)任何處理的原始圖像，圖4(c)為線性差值獲得的結(jié)果，圖4(d)、(e)、(f)為本發(fā)明方法一次迭代、二次迭代、三次迭代后的結(jié)果圖，圖(g)、(h)、(i)為patch方法一次迭代、二次迭代、三次迭代后的結(jié)果圖。具體實施方式實施例1。一種圓軌道下不完全ct投影數(shù)據(jù)的恢復方法，通過如下步驟進行：第一步，獲得壞像素點位置信息，將實際采集的投影數(shù)據(jù)轉(zhuǎn)換成對應的弦圖，即：弦圖的x軸為投影的像素點位置，y軸為該像素位置對應的投影角度位置；對某個探測器單元存在損壞的條件下，則在弦圖中對應的像素位置上會出現(xiàn)一條豎線，根據(jù)此理論在弦圖中獲得壞像素點的位置信息，并記錄壞像素點位置。第二步，對每個投影中的壞像素點使用水平一維三次樣條插值方法進行線性插值，以恢復壞像素點的低頻信息，經(jīng)過此步驟得到每個角度下的初始化投影。第三步，令n＝1，進入第四步，第四步，令m＝n，m為當前迭代次數(shù)；將每個角度下的初始化投影進行二維傅里葉變換，獲得每個角度下的初始化投影在頻率域中的頻域數(shù)據(jù)；第五步，分別將前后兩個相鄰投影的頻率域形式帶入到一致性條件對應的公式中計算投影結(jié)果；具體是使用相同的權(quán)重因子ω對前后兩個投影結(jié)果進行加權(quán)，得到加權(quán)更新后的投影結(jié)果。第五步中使用的一致性條件，形式如下：其中，其中ρ是射線源到旋轉(zhuǎn)中心的距離，d是旋轉(zhuǎn)中心到探測器的距離，(α1,α2)是平板探測器坐標，θ是方位角，如圖1所示。第五步，可具體使用一致性條件即公式(28)對每一個初始化投影的頻域數(shù)據(jù)進行恢復，公式(28)的具體形式如下：u*(θ+dθ)＝u*(θ)+dθ·u,k2≠0……公式(28)；其中，u*(θ+dθ)和u*(θ)分別表示θ+dθ和θ角度下初始化投影的頻域數(shù)據(jù)，u*(θ)為已知量，u*(θ+dθ)是要求的未知量，dθ代表相鄰兩個角度的間隔，k1、k2是平板探測器的頻率域坐標，i代表復數(shù)的虛部；公式(28)的具體使用方法如下：首先將u*(θ+dθ)代入到公式(28)中計算u*(θ)：u*(θ)＝u*(θ+dθ)-dθ·u；然后用u*(θ-dθ)代替公式(28)中的u*(θ)，u*(θ)代替公式(28)中的u*(θ+dθ)計算得到：u*(θ)＝u*(θ-dθ)+dθ·u；由于選取前后兩個相鄰投影頻域數(shù)據(jù)的任一個都可以實現(xiàn)投影數(shù)據(jù)恢復，為了平衡恢復效果，使用權(quán)重因子ω，ω∈[0,1]，對上述兩個公式的計算結(jié)果進行加權(quán)平均，得到加權(quán)后的結(jié)果ur*(θ)，即:ur*(θ)＝ω(u*(θ-dθ)+dθ·u)+(1-ω)(u*(θ+dθ)-dθ·u),k2≠0；由于公式(28)只有在k2≠0時才有效，所以對頻域數(shù)據(jù)u*(θ)中k2＝0的部分，直接使用θ角度下初始化投影頻域數(shù)據(jù)中的相應部分來補償；第六步，具體是對ur*(θ)進行傅里葉反變換，得到第m次迭代終值ur(θ)；ur(θ)＝f-1(ur*(θ))……(1)；其中公式(1)中ur(θ)表示θ角度下加權(quán)后的頻域數(shù)據(jù)ur*(θ)的反傅里葉變換結(jié)果，f-1代表傅里葉反變換。第七步，具體是判斷當前迭代次數(shù)m是否大于等于迭代終止次數(shù)s，如果是，則進入第九步；否則進入第八步。第八步，具體是將第m次迭代終值ur(θ)作為初始化投影，令n＝n+1，進入第三步。第九步，具體是以第m次迭代終值ur(θ)為最終恢復的數(shù)據(jù)結(jié)果。其中本方法中應用的公式(28)推導過程如下：由于本發(fā)明提出的一致性條件是基于約翰方程，所以先對約翰方程進行簡單介紹和分析。fjohn于1938年在文獻《theultrahyperbolicdifferentialequationwithfourindependentvariables》中提出一個超雙曲線二階偏微分方程，作為線積分的一致性條件，具體形式如下：u是函數(shù)f穿過ε和η兩點的線積分：u(ε；η)＝∫f(ε+t(η-ε))dt(3)；以下為圓軌道下基于約翰方程的一致性條件推導步驟：首先采用圓軌道cbct坐標參數(shù)對ε和η進行參數(shù)化：其中ρ是射線源到旋轉(zhuǎn)中心的距離，d是旋轉(zhuǎn)中心到探測器的距離，(α1,α2)是平板探測器坐標，θ是方位角。分別求ε和η的每個分量求偏導：為書寫方便，定義如下算子：把公式(6)-(11)代入到公式(2)中得到：由公式(13)-(16)對h1，h2,的定義，代入到公式(17)-(19)中，得到：從上述公式中，得出：對上面公式(21)移項合并同類項：求解公式(23)，最終推導結(jié)果如下：公式(24)即為我們推導的圓軌道下基于約翰方程的一致性條件。對公式(24)進行傅里葉變換，得到對應的頻率域的表達形式如下：在上述公式中，uθ*表示u對θ求導后的傅里葉變換，(k1,k2)是(α1,α2)的對應頻率域形式。為了方便表達和應用，將公式(25)右邊式子用u代替：考慮到求導公式的一般定義，可以將公式(25)改下成如下形式：即：u*(θ+dθ)＝u*(θ)+dθ·u,k2≠0(28)；公式(28)是本發(fā)明改進后的一致性條件，本發(fā)明使用公式(28)進行投影數(shù)據(jù)的恢復。傳統(tǒng)的空間插值技術(shù)使用同一個角度的投影的其他完整像素點數(shù)據(jù)來恢復丟失的數(shù)據(jù)，所使用重建信息只局限于丟失數(shù)據(jù)所在角度的投影，所以插值方法只能恢復丟失數(shù)據(jù)中的低頻信息，對丟失的高頻信息基本沒有恢復。本發(fā)明方法應用的一致性條件，即不僅基于同一個投影中的數(shù)據(jù)來恢復丟失數(shù)據(jù)，還利用相鄰角度投影中的對應位置數(shù)據(jù)信息。因此本發(fā)采用一致性條件可以很好地恢復丟失數(shù)據(jù)中的高頻信息，能基本消除使用插值技術(shù)獲得的重建圖像中的條形偽影，顯著提高圖像質(zhì)量。和傳統(tǒng)的空間插值方法相比，本發(fā)明的方法能夠很好地恢復丟失信息中的高頻成分，進而大幅度減少圖像中的條形偽影，顯著提高圖像質(zhì)量。在圓軌道下使用patch的方法需要先使用探測器的縱坐標來近似代替z坐標，而本方法采用的一致性條件是基于圓軌道cbct幾何結(jié)構(gòu)推導得到的，能夠避免因為近似z坐標而導致的精度下降，并且patch的方法中z軸在不停變化，使得此方法具體應用困難。與之相反，本發(fā)明方法只是一個簡單的一階求導公式，實現(xiàn)簡單，具有可應用性。實施例2。本實例分別使用傳統(tǒng)空間插值方法、patch的方法(
背景技術(shù)：
中有詳細說明)、以及本方法對一組不完整的投影數(shù)據(jù)進行恢復，不完整的投影數(shù)據(jù)由一臺平板探測器具有壞像素點的cbct采集得到，其ct掃描參數(shù)具體如下：射線源到旋轉(zhuǎn)中心的距離ρ為500mm，旋轉(zhuǎn)中心到探測器的距離d為500mm，平板探測器大小為850×200mm，像素點尺寸為1×1mm，圍繞物體旋轉(zhuǎn)360度采集1080個投影。本發(fā)明技術(shù)的方法具體實施方式如下：第一步，獲得壞像素點位置信息。將實際采集的投影數(shù)據(jù)轉(zhuǎn)換成對應的弦圖，如圖2所示。圖中黑色豎線所對應的橫坐標即為壞像素點位置，通過此方式可以找到所有壞像素點的位置信息，并記錄。第二步，對每一個投影中的壞像素點使用水平一維三次樣條插值方法線性插值，并將插值后的結(jié)果作為該壞像素點的初始值。第三步，將插值后的投影以及前后兩個相鄰投影進行傅里葉變換，得到對應的頻率域形式，以方便應用本發(fā)明推導出的一致性條件。第四步，將上一步得到的前后兩個相鄰投影的頻率域形式代入到公式(28)中，恢復中間投影的丟失信息(尤其是高頻成分信息)，并使用同一個權(quán)重因子ω對兩個結(jié)果進行加權(quán)，由于沒有更多的約束條件表明前后兩個相鄰投影哪個對中間投影的貢獻更大，所以本實例令ω＝0.5。然后用計算加權(quán)后的結(jié)果修正該投影初始值。第五步，對修正后的值進行傅里葉反變換，得到第一次迭代最終值。第六步，將上一次迭代的最終結(jié)果作為初值，重復上面第三步到第五步，本實例共迭代三次。具體步驟流程圖如圖3所示。仿真結(jié)果如圖4所示，圖4(a)為理想圖像，圖4(b)為未經(jīng)任何處理的原始圖像，圖4(c)為線性差值獲得的結(jié)果，圖4(d)、(e)、(f)為本發(fā)明方法一次迭代、二次迭代、三次迭代后的結(jié)果圖，圖(g)、(h)、(i)為patch方法一次迭代、二次迭代、三次迭代后的結(jié)果圖。使用絕對誤差(absoluteerror)對重建圖像進行分析:其中n為像素點個數(shù)，ri是重建圖像第i個點對應的灰度值，ti是理想圖像第i個點對應的灰度值。結(jié)果如下：方法第101層的mae線性空間插值1.3368×10+6patch方法一次迭代3.0641×10+5patch方法二次迭代3.1477×10+5patch方法三次迭代2.6800×10+5本方法一次迭代2.9724×10+5本方法二次迭代2.9658×10+5本方法三次迭代2.6187×10+5從重建圖像效果來看，使用空間插值方法獲得的重建圖像條形偽影嚴重，嚴重影響了圖像診斷信息的提取，而patch方法和本方法則可以顯著地改善圖像質(zhì)量。使用絕對誤差對重建圖像進行定量分析可以發(fā)現(xiàn)：和patch方法相比，本發(fā)明方法獲得的圖像質(zhì)量更接近理想圖像質(zhì)量，每一次的迭代效果都優(yōu)于patch方法，而且后面的迭代效果優(yōu)于前一次的迭代效果。最后應當說明的是，以上實施例僅用以說明本發(fā)明的技術(shù)方案而非對本發(fā)明保護范圍的限制，盡管參照較佳實施例對本發(fā)明作了詳細說明，本領域的普通技術(shù)人員應當理解，可以對本發(fā)明的技術(shù)方案進行修改或者等同替換，而不脫離本發(fā)明技術(shù)方案的實質(zhì)和范圍。當前第1頁12