專利名稱：一種基于球諧函數的磁共振圖像梯度變形校正方法
技術領域：
本發明涉及一種基于球諧函數的磁共振圖像梯度變形校正方法。
背景技術：
梯度場的非線性會引起圖像的幾何變形，在不同的磁共振成像設備上以及同樣的磁共振成像設備而不同的視野(FOV)內圖像的幾何變形是不同的。針對梯度場非線性的圖像校正對提高磁共振圖像導航手術的精度和多磁共振之間圖像研究的可行性研究是必要的。磁共振圖像梯度變形校正的基本過程是先通過某種方法求出磁共振成像設備成像空間內有限個控制點的位置偏移量，接下來計算待校正圖像的每一個象素在磁共振成像設備成像空間內的坐標，通過插值方法根據每一個象素點的坐標計算該象素的空間偏移量，根據磁共振成像設備坐標系與圖像坐標系的轉換關系，將該象素的空間偏移量換算為圖像坐標系內的偏移量，將該偏移量補償到該象素的圖像坐標，即完成了磁共振圖像梯度變形校正過程。計算磁共振成像空間內控制點的位置偏移量的方法可以分為兩種第一種是使用三維立體水模，三維立體水模提供了三維空間內的控制點，通過對立體水模的磁共振圖像內控制點的分析，可以得到有限空間內控制點上的圖像變形量，根據磁共振成像設備坐標系與圖像坐標系的轉換關系，將該象素的圖像坐標系內的偏移量換算為空間偏移量。另一種方法是應用磁場的描述函數-球諧函數計算成像空間內有限個控制點的偏移量。由于第一種方法依賴于三維水模，實施成本遠高于第二種方法。而且需要對三維水模多次成像和分析才能獲得控制點的位置偏移量，所以實施的效率很低。第二種方法在應用球諧函數進行計算之前，首先要給出球諧函數內的相關校正參數。哈佛大學的研究人員Jorge Jovicich，Silvester Czanner提供了第二種校正方法以及程序，論文題目為Reliability in Multi-SiteStructural MRI StudiesEffects of Gradient Non-linearity Correction onPhantom and Human Data，NeuroImage(in press)(多點結構性磁共振成像可靠性研究梯度非線性校正在水模和人體圖像中的效果，《神經圖像》中文譯文)。該方法是根據梯度線圈的設計參數計算用于圖像校正的校正參數，但梯度線圈的設計參數只是圖像變形的一個主要因素，并沒有考慮其它因素對圖像變形所產生的影響，所以實際的校正參數與設計參數仍存在一定差異。

發明內容
針對上述問題，本發明的主要目的在于提供一種針對特定的梯度線圈在圖像空間內取有限的點精確計算校正參數，從而獲取真實的磁場梯度參數的基于球諧函數的磁共振圖像梯度變形校正方法。
為達到上述目的，本發明采取以下技術方案一種基于球諧函數的磁共振圖像梯度變形校正方法，其特征在于它包括計算空間偏移量，圖像偏移量補償校正和搜索校正參數三個步驟，其是根據磁感應強度的球諧函數，計算三維空間內的圖像上每一個象素點的三維偏移量，并將該偏移量補償到該象素點的圖像坐標中。
所述計算空間偏移量的步驟包括(a)使用公式(1)計算在球坐標系內磁共振系統的磁場強度Bz(r，θ，φ)Br(n，m)(r，θ，φ)＝rn[av(n，m)cos(mφ)+bv(n，m)sin(mφ)]×P(n，m)(cosθ) (1)(b)得到磁場的梯度函數Gv(r)≡∂Bzv(r)∂v≡∂Bzv(r)L∂v+∂Bzv(r)N∂v≡GvL+Gv(r)N---(2)]]>Bzv(r)N是可以由球諧函數展開的高階項計算的非線性部分，所以 定義 其中 根據磁共振成像設備的成像空間選擇計算區間，并在該區間內均勻分布控制點Vci(x，y，z)，根據公式(4)計算每一個控制點的偏移量，并保存在控制點偏移量文件中。
所述圖像偏移量補償校正的步驟包括(a)將象素p(u，v)的圖像坐標(u，v)換算為成像空間內坐標(x，y，z)；(b)根據控制點偏移量文件中所保存的與該象素相鄰的8個控制點Vc(1-8)的偏移量ηc(1-8)使用插值算法計算該象素的空間偏移量ηp；(c)將該象素的空間偏移量換算為圖像坐標系內的偏移量σp；(d)將偏移量補償到該象素的圖像坐標中。
所述搜索校正參數的步驟包括a)梯度場線性標定在梯度磁場強度為線性函數的視野內，通過掃描已知尺寸的水模，求出物理尺寸與圖像坐標尺寸的比例因子；b)設置初始的系統參數Cf設定X方向ax(5，1)、ax(3，1)初始參數Cxf；c)計算空間位置偏移量ηx將設定的初始參數Cxf帶入公式 中，計算出空間位置偏移量ηx，并存入偏移量查找表；d)將梯度校正水模放置在磁場的中心附近，并保證水模圓柱的行的方向與X方向一致，列方向與Y或Z方向一致；e)設置磁共振儀器的掃描參數，在步驟d)的放置狀態下，在x，y平面內設置掃描的層，圖像中心位于磁場中心，視野為希望被校正的范圍，并得到梯度校正水模磁共振原始圖像Mo；f)使用圖像偏移量補償校正模塊對原始圖像校正，得到校正后的水模圖像Mc；g)取包含x軸上的標志點p1，2，3…的圖像進行處理，采用重心法計算每一個標志點的圖像坐標Mi；h)根據步驟a)求出的比例因子計算標志點p1，2，3…對應的世界坐標Xi；i)計算標志點p1，2，3…的位置誤差平方和S；j)判斷S是否滿足要求或迭代次數n是否超過設定值，如果“是”則退出迭代，以當前系統參數值為程序輸出；如果“否”則改變初值繼續進行搜索。
所述步驟e)中磁共振儀器的掃描參數包括視野、掃描方向，圖像中心的位置。
所述校正水模按照磁場梯度變化的規律采用不等間距的標志點分布和尺寸不同的標志點，標志點間距和尺寸的計算方法如下根據公式(3)以及ax(n，m)、bx(n，m)的估計值ax(3，1)＝10-4，ax(5，1)＝10-7可以寫出v(x，y＝0，z＝0)處的梯度場磁感應強度Bz的5階展開BzX(V)=x-32*10-4*x3+316*10-6*x5]]>v(x，y＝0，z＝0)處的梯度為GX(v)=1-92*10-4*x2+332*10-5*x4]]>假設磁場中心的標志點間距為l0，標志點的直徑為d0。那么的v(x，y＝0，z＝0)點處的標志點間距為l=l0G;]]>v(x，y＝0，z＝0)點處的標志點直徑為d=d0GX(v)]]>本發明由于采取以上技術方案，其具有以下優點1、本發明根據實際圖像反求磁共振校正參數，克服了梯度線圈設計參數與真實的磁場梯度參數之間存在差異的缺點。2、本發明在反求磁場梯度參數的過程中，采取依次反求三個方向的磁場梯度參數的方法，提高了算法的收斂速度。3、本發明采用校正水模非均勻標志點分布方案，有效地提高了校正精度。本發明方法可以廣泛用于磁共振圖像引導下的手術治療過程中，對病灶圖像形成的梯度變形的校正中，使經本發明校正后的圖像質量更加真實有效。


圖1是本發明構成原理示意2是本發明水模結構示意3是圖2的側視剖視示意4是本發明X軸上的磁場梯度在±30cm內的變化規律示意5是本發明等間距標志點的圖像變化示意6是按照標志點的間距和直徑計算公式設計的水模標志點示意7是本發明的流程8是放置水模的一實施方式示意9是梯度校正水模磁共振原始圖像Mo圖10是本發明實施例校正后的水模圖像圖11是本發明包含x軸上的標志點p1，2，3…示意圖具體實施方式
下面結合附圖和實施例對本發明進行詳細的描述。
本發明所提供的基于球諧函數的磁共振圖像梯度變形校正方法，是根據磁感應強度的球諧函數，計算三維空間內的圖像上每一個象素點的三維偏移量，并將該偏移量補償到該象素點的圖像坐標中。
如圖1所示，本發明包括計算空間偏移量、圖像偏移量補償校正和搜索校正參數三部份。
1、計算空間偏移量的步驟如下(a)磁共振系統的磁場強度Bz(r，θ，φ)在球坐標系內可以使用公式1表示。
Br(n，m)(r，θ，φ)＝rn[av(n，m)cos(mφ)+bv(n，m)sin(mφ)]×P(n，m)(cosθ)(1)其中Br(n，m)(r，θ，φ)是Bz的球諧函數展開后第n階m級項。av(n，m)、bv(n，m)是常數，av(n，m)、bv(n，m)是v方向n階m級展開項的系數，是磁場非線性梯度的固有特性。r為所求點到球心的距離。P(n，m)(cosθ)為勒讓德多項式。有限項Br(n，m)(r，θ，φ)之和就可以近似給出磁場強度Bz(r，θ，φ)。
(b)在建立了Bz(r，θ，φ)的表達式之后，就可以得到磁場的梯度函數Gv(r)≡∂Bzv(r)∂v≡∂Bzv(r)L∂v+∂Bzv(r)N∂v≡GvL+Gv(r)N---(2)]]>其中v表示笛卡兒坐標系內的(x，y或z)方向，Bzv(r)是梯度場的總強度，Bzv(r)L是梯度場線性部分，即球諧函數展開的一階項，Bzv(r)N是可以由球諧函數展開的高階項計算的非線性部分，所以 定義 其中
公式(4)描述了笛卡兒坐標系內v(x，y或z)方向上梯度的非線性程度，可以用于計算笛卡兒空間內點V(x，y，z)的坐標偏移量，可以根據磁共振成像設備的成像空間選擇合適大小的計算區間(一般能夠覆蓋MRI在掃描過程中所使用的所有成像空間的最小體積，例如掃描人體最大的體積也就是0.6m×0.6m×0.6m)，并在該區間內均勻分布控制點Vci(x，y，z)，根據公式(4)計算每一個控制點的偏移量，并保存在控制點偏移量文件中。
2、偏移量補償校正的步驟如下對待校正圖像進行校正由圖像偏移量補償校正實現，是在圖像坐標系中通過改變每一個象素的圖像坐標完成的，所以針對每一個象素的圖像校正的過程，即(a)首先將象素p(u，v)的圖像坐標(u，v)換算為成像空間內坐標(x，y，z)；(b)再根據控制點偏移量文件中所保存的與該象素相鄰的8個控制點Vc(1-8)的偏移量ηc(1-8)使用插值算法計算該象素的空間偏移量ηp；(c)再次將該象素的空間偏移量換算為圖像坐標系內的偏移量σp；(d)最后將偏移量補償到該象素的圖像坐標。具體公式描述如下假設在磁共振成像空間內待校正圖像的行對應的向量為vr，列對應的向量為vc。第一個象素點為P0(u0，v0)，其磁共振空間坐標為V0，圖像的分辨率為Res，圖像的視野為Fov。則任意象素Pi(u，v)的磁共振空間坐標為Vi。
Vi(xi,yi,zi)=V0+((u-u0)*vr+(v-v0)*vc)*FovRes---(5)]]>假設與該象素相鄰的8個控制點的偏移量分別為η(x0，y0，z0)，η(x0，y0，z0+1)，η(x0，y0+1，z0)，η(x0+1，y0，z0)，η(x0，y0+1，z0+1)，η(x0+1，y0+1，z0)，η(x0+1，y0，z0+1)，η(x0+1，y0+1，z0+1)。
則Vi(xi，yi，zi)的偏移量可以表示為ηi(xi,yi,zi)=η(x0,y0,z0)*(x0+1-xi)*(y0+1-yi)*(z0+1-zi)]]>+η(x0+1,y0,z0)*(xi-x0)·*(y0+1-yi)*(z0+1-zi)]]>+η(x0,y0+1,z0)*(x0+1-xi)*(yi-y0)*(z0+1-zi)]]>+η(x0,y0,z0+1)*(x0+1-xi)*(y0+1-yi)*(zi-z0)]]>+η(x0+1,y0,z0+1)*(xi-x0)*(y0+1-yi)*(zi-z0)]]>
+η(x0,y0+1,z0+1)*(x0+1-xi)*(yi-y0)*(zi-z0)]]>+η(x0+1,y0+1,z0)*(xi-x0)*(yi-y0)*(z0+1-zi)]]>+η(x0+1,y0+1,z0+1)*(xi-x0)*(yi-y0)*(zi-z0)---(6)]]>將該象素的空間偏移量ηi(xi，yi，zi)換算為圖像坐標系內的偏移量σp(u，v)。
σp(u，v)＝(Δu，Δv)Δu=ηi(xi,yi,zi)·vr*ResFov]]>Δv=ηi(xi,yi,zi)·vc*ResFov---(7)]]>將偏移量補償到該象素的圖像坐標，該象素新的圖像坐標(u′，v′)。
u′＝u+Δu；v′＝v+Δv；(8)對待校正圖像內的每一個象素進行以上操作，就可以完成整幅圖像的校正。
3、搜索校正參數搜索校正參數模塊通過多次分析梯度校正水模的磁共振圖像獲得校正參數，下面介紹水模的設計方法和搜索校正參數的步驟，出于精度和計算量的考慮梯度場的計算均以5階展開。
(1)梯度校正水模如圖2、圖3所示，梯度校正水模是用于梯度校正的工具，它既用于搜索系統參數，也用于測量校正后的誤差。梯度校正水模為一個方形的盒體1，在盒體1內設置按照不等間距方陣的方式排列成方陣的圓柱2，圓柱2內充滿硫酸銅溶液，可以在磁共振成像設備中成像。梯度校正水模中的圓柱2是按照公式(2)表達的梯度分布規律進行排列的，圓柱2作為水模內的標志點，可以更精確地計算出梯度分布的參數值。
如圖4所示，顯示了圓柱2在X軸上的磁場梯度在±30cm內的變化規律，由于磁場梯度的非線性，以X方向標志點為例，如果標志點為等間距的，則隨著X坐標越大，相鄰標志點間的圖像X坐標差值先變大，之后急劇變小(如圖5所示)。所以為了保證標志點的圖像坐標計算精度，本發明按照磁場梯度變化的規律采用不等間距的標志點分布和尺寸不同的標志點。下面以X方向為例給出標志點間距和尺寸的計算方法根據公式(3)以及ax(n，m)、bx(n，m)的估計值ax(3，1)＝10-4，ax(5，1)＝10-7，可以寫出v(x，y＝0，z＝0)處的梯度場磁感應強度Bz的5階展開
BzX(v)=x-32*10-4*x3+316*10-6*x5]]>v(x，y＝0，z＝0)處的梯度為GX(v)=1-92*10-4*x2+332*10-5*x4]]>假設磁場中心的標志點間距為l0，標志點的直徑為d0。那么的v(x，y＝0，z＝0)點處的標志點間距為l=l0G;]]>v(x，y＝0，z＝0)點處的標志點直徑為d=d0GX(v)]]>同樣可以根據公式(3)以及az(n，m)、bz(n，m)的估計值得出梯度場磁感應梯度在z軸方向上的表達式GZ。進而根據GZ計算標志點在平面水模另一個方向的點距以及標志點尺寸。一般來講磁共振的梯度場x，y方向的設計參數基本相同，所以x、y方向僅計算一個即可，一個典型的梯度校正用平面水模可以使用x，z方向的分布數據制作。
按照標志點的間距和直徑計算公式重新設計本發明水模(如圖6所示)。
(2)搜索校正參數的步驟如下校正參數搜索模塊需要分三次分別搜索X、Y和Z方向的校正參數，下面僅以X方向校正參數搜索為例介紹實施步驟(如圖7所示)a)梯度場線性標定假設在5cm的視野內，梯度磁場強度為線性函數，通過掃描已知尺寸的水模，求出物理尺寸與圖像坐標的比例因子；b)設置初始的系統參數Cf設定X方向ax(5，1)、ax(3，1)初始參數Cxfc)計算空間位置偏移量ηx將設定的初始參數Cxf帶入公式 中，計算出空間位置偏移量ηx，并存入偏移量查找表。
d)將梯度校正水模放置在磁場的中心附近，并保證水模圓柱的行的方向與X方向一致，列方向與Y或Z方向一致(如圖8所示)。
e)設置磁共振儀器的掃描參數，包括視野、掃描方向，圖像中心的位置，設置掃描的層(slice)在x，y平面內，圖像中心位于磁場中心，視野為希望被校正的范圍，并可得到梯度校正水模磁共振原始圖像Mo(如圖9所示)。
f)使用圖像偏移量補償校正模塊對原始圖像校正，得到校正后的水模圖像Mc(如圖10所示)。
g)取包含X軸上的標志點p1，2，3…的圖像(如圖11所示)進行處理，采用重心法計算每一個標志點的圖像坐標Mi。
h)根據步驟a)求出的比例因子計算標志點p1，2，3…對應的世界坐標Xi。
i)計算標志點p1，2，3…的位置誤差平方和S。
j)判斷S是否滿足要求或迭代次數n是否超過設定值，如果“是”則退出迭代，以當前系統參數值為程序輸出；如果“否”則改變初值繼續進行搜索。
采用上述校正方法，在梯度線圈視野為5cm內，梯度能夠保證線性度誤差在4‰左右。
權利要求
1.一種基于球諧函數的磁共振圖像梯度變形校正方法，其特征在于它包括計算空間偏移量，圖像偏移量補償校正和搜索校正參數三個步驟，其是根據磁感應強度的球諧函數，計算三維空間內的圖像上每一個象素點的三維偏移量，并將該偏移量補償到該象素點的圖像坐標中。
2.如權利要求1所述的一種基于球諧函數的磁共振圖像梯度變形校正方法，其特征在于所述計算空間偏移量的步驟包括(a)使用公式(1)計算在球坐標系內磁共振系統的磁場強度Bz(r，θ，φ)Br(n，m)(r，θ，φ)＝rn[av(n，m)cos(mφ)+bv(n，m)sin(mφ)]×P(n，m)(cosθ) (1)(b)得到磁場的梯度函數Gv(r)≡∂Bzv(r)∂v≡∂Bzv(r)L∂v+∂Bzv(r)N∂v≡GvL+Gv(r)N---(2)]]>Bzv(r)N是可以由球諧函數展開的高階項計算的非線性部分，所以 定義 其中 根據磁共振成像設備的成像空間選擇計算區間，并在該區間內均勻分布控制點Vci(x，y，z)，根據公式(4)計算每一個控制點的偏移量，并保存在控制點偏移量文件中。
3.如權利要求1所述的一種基于球諧函數的磁共振圖像梯度變形校正方法，其特征在于所述圖像偏移量補償校正的步驟包括(a)將象素p(u，v)的圖像坐標(u，v)換算為成像空間內坐標(x，y，z)；(b)根據控制點偏移量文件中所保存的與該象素相鄰的8個控制點Vc(1-8)的偏移量ηc(1-8)使用插值算法計算該象素的空間偏移量ηp；(c)將該象素的空間偏移量換算為圖像坐標系內的偏移量σp；(d)將偏移量補償到該象素的圖像坐標中。
4.如權利要求2所述的一種基于球諧函數的磁共振圖像梯度變形校正方法，其特征在于所述圖像偏移量補償校正的步驟包括(a)將象素p(u，v)的圖像坐標(u，v)換算為成像空間內坐標(x，y，z)；(b)根據控制點偏移量文件中所保存的與該象素相鄰的8個控制點Vc(1-8)的偏移量ηc(1-8)使用插值算法計算該象素的空間偏移量ηp；(c)將該象素的空間偏移量換算為圖像坐標系內的偏移量σp；(d)將偏移量補償到該象素的圖像坐標中。
5.如權利要求1或2或3或4所述的一種基于球諧函數的磁共振圖像梯度變形校正方法，其特征在于所述搜索校正參數的步驟包括a)梯度場線性標定在梯度磁場強度為線性函數的視野內，通過掃描已知尺寸的水模，求出物理尺寸與圖像坐標尺寸的比例因子；b)設置初始的系統參數Cf設定X方向ax(5，1)、ax(3，1)初始參數Cxf；c)計算空間位置偏移量ηx將設定的初始參數Cxf帶入公式 中，計算出空間位置偏移量ηx，并存入偏移量查找表；d)將梯度校正水模放置在磁場的中心附近，并保證水模圓柱的行的方向與X方向一致，列方向與Y或Z方向一致；e)設置磁共振儀器的掃描參數，在步驟d)的放置狀態下，在x，y平面內設置掃描的層，圖像中心位于磁場中心，視野為希望被校正的范圍，并得到梯度校正水模磁共振原始圖像Mo；f)使用圖像偏移量補償校正模塊對原始圖像校正，得到校正后的水模圖像Mc；g)取包含x軸上的標志點p1，2，3…的圖像進行處理，采用重心法計算每一個標志點的圖像坐標Mi；h)根據步驟a)求出的比例因子計算標志點p1，2，3…對應的世界坐標Xi；i)計算標志點p1，2，3…的位置誤差平方和S；j)判斷S是否滿足要求或迭代次數n是否超過設定值，如果“是”則退出迭代，以當前系統參數值為程序輸出；如果“否”則改變初值繼續進行搜索。
6.如權利要求5所述的一種基于球諧函數的磁共振圖像梯度變形校正方法，其特征在于所述步驟e)中磁共振儀器的掃描參數包括視野、掃描方向，圖像中心的位置。
7.如權利要求5所述的一種基于球諧函數的磁共振圖像梯度變形校正方法，其特征在于所述校正水模按照磁場梯度變化的規律采用不等間距的標志點分布和尺寸不同的標志點，標志點間距和尺寸的計算方法如下根據公式(3)以及ax(n，m)、bx(n，m)的估計值ax(3，1)＝10-4，ax(5，1)＝10-7可以寫出v(x，y＝0，z＝0)處的梯度場磁感應強度Bz的5階展開BzX(v)=x+32*10-4*x3+310*10-6*x5]]>v(x，y＝0，z＝0)處的梯度為GX(v)=1-92*10-4*x2+332*10-5*x4]]>假設磁場中心的標志點間距為l0，標志點的直徑為d0。那么的v(x，y＝0，z＝0)點處的標志點間距為l=l0G;]]>v(x，y＝0，z＝0)點處的標志點直徑為d=d0GX(v).]]>
全文摘要
本發明涉及一種基于球諧函數的磁共振圖像梯度變形校正方法，其特征在于它包括計算空間偏移量，圖像偏移量補償校正和搜索校正參數三個步驟，其是根據磁感應強度的球諧函數，計算三維空間內的圖像上每一個象素點的三維偏移量，并將該偏移量補償到該象素點的圖像坐標中。本發明根據實際圖像反求磁共振校正參數，克服了梯度線圈設計參數與真實的磁場梯度參數之間存在差異的缺點；同時在反求磁場梯度參數的過程中，采取依次反求三個方向的磁場梯度參數的方法，提高了算法的收斂速度；另外本發明采用校正水模非均勻標志點分布方案，有效地提高了校正精度。本發明方法可以廣泛用于磁共振圖像引導下的手術治療過程中，對病灶圖像形成的梯度變形的校正中，使經本發明校正后的圖像質量更加真實有效。
文檔編號G01R33/58GK101046506SQ200710064899
公開日2007年10月3日 申請日期2007年3月28日 優先權日2007年3月28日
發明者趙磊, 代亮 申請人:新奧博為技術有限公司