一種基于線性關系的熒光分子斷層成像重建方法

專利名稱：：一種基于線性關系的熒光分子斷層成像重建方法
技術領域：
：本發明屬于在生物醫學中圖像重建的
技術領域：
，涉及一種三維圖像重建方法。技術背景分子成像是近年來生物醫學領域的一個熱門方向。分子成像，廣義上可以定義為在分子和細胞水平對活體狀態下的生物過程進行的定性和定量的測量[R.Weissleder，Radiology,2001]。熒光分子斷層成像(FMT)是分子成像中的一個重要方面。它被廣泛的應用在生物，醫學研究領域，使得熒光分子成像已經從一種定性的用作切片檢査的工具變成了一種能夠用于準確的，3維成像的工具。熒光分子斷層成像的關鍵在于獲得多個角度的投影數據用于重建。同CT等基于高能射線的設備有所不同，光(400nm900nm)在生物組織中的傳播不是遵循直線傳播。這個光傳播路徑可以近似用漫射方程來描述。利用基于漫射方程模型的重建方法，從多個角度投影數據就可以斷層地重建出熒光探針的濃度分布。R.B.Schulz等人[R.B.Schulz，正EETransactionsonMedicalImaging,2004]提出了一個基于解析解的重建方法，能夠得到表面檢測到熒光強度同內部熒光體濃度分布之間的線性關系。這種方法假設成像物體的光學參數是均勻的。然而，實際成像的實驗用鼠的光學參數是非均勻的。因而，光學參數估計的不準確性會導致重建結果的很大失真。S.Bjoem等人[S.Bjoem，AdvancesinMedicalEngineering,2007]對光學參數估計不準確造成失真做了研究。白凈教授領導的醫學成像實驗室在2007年發表了基于有限元方法的重建方法[X.Song，OpticsExpress,2007]，能夠處理非均勻的光學參數，同時在方法中計算得到表面檢測到的熒光強度同未知的熒光探針分布之間的線性關系，進而進行重建。這些方法具有準確，在迭代過程中不需要反復求解梯度以及二屆梯度的特點，因而計算速度以5-10分鐘計。然而，在計算線性關系的過程中，有針對大規模稀疏矩陣求逆的操作，限制了方法的計算速度地進一步提高
發明內容本發明的目的是為熒光分子斷層成像系統提供一種重建方法，其特點是能夠快速準確地建立實驗用鼠表面檢測到的熒光強度同未知的熒光探針分布之間的線性關系，實現重建。在這個過程中，方法能夠處理不規則的成像物體以及非均勻的光學參數。方法中在其它成像模式，如CT，得到的分割好的實驗用鼠圖譜的基礎上，利用二值圖像輪廓線提取以及數值擬合的方法得到實驗用鼠的不規則外表面輪廓；在實驗用鼠圖譜上，通過査表的方法確定實驗用鼠體內的非均勻光學參數，包括吸收和散射參數；針對每個激發光源-熒光檢測對，利用一個向量-稀疏矩陣乘積的方法得到這個檢測點處的熒光強度同未知的熒光探針分布之間的線性關系。本發明的特征在于，所述方法是用一臺計算機依次按以下步驟實現的-.步驟(l).把用CT成像設備分割完畢的生物實驗用鼠的分層的圖譜輸入所述計算機中，該計算機通過下述圖像處理方法得到外表面輪廓，并把所表述表面輪廓建模成幾何體步驟(l.l).對所述圖譜的每--層圖像/，用一個相同的閾值*=0.5進行二值化大于所選閾值Ar者取l，小于或等于該閾值Ar者取O，把該圖譜中的組織和背景依次用所述的"1"和"0"分開，得到每一層圖像的二值圖BI:<formula>formulaseeoriginaldocumentpage8</formula>其中，row代表該圖像/的行，co/代表該圖像/的列，步驟(1.2).用二值圖輪廓線提取的方法得到所述的二值圖S/的外輪廓線S:對該二值圖W上的每一個像素點5/[raw，co/]執行以下操作若5/[row,co/]為1,而且其相鄰的8個像素點至少有一個為0，則該像素點萬/[raw，co/]在所述外輪廓線S上，若S/[row,co/]為0，而且其相鄰的8個像素點至少有一個為1，則該像素點B/[row,co/]在所述外輪廓線S上，步驟(1.3).把步驟(1.2)得到的所述外輪廓線S從像素坐標對應到笛卡爾坐標，則該外輪廓線S上的任意一個點[raw,co/]對應到所述笛卡爾坐標[je，W上為<formula>formulaseeoriginaldocumentpage8</formula>其中，點[row。，co/。]為該輪廓線的中心，A是每個像素點的幾何尺寸，用cw表示，步驟(1.4).把步驟(1.3)得到的該笛卡爾坐標[x，y]對應到極坐標[cr,i]上，《是角度，i是極坐標的半徑，對所述外輪廓線S上的所有點，把每個點上的所述半徑i用角度a按下式做多項式擬合<formula>formulaseeoriginaldocumentpage9</formula>c為擬合次數，c=0，l,'.'12，A是擬合得到的多項式系數，用{八.}表示，"二l,2,…&m^，S—m^是所述輪廓線S上的所有點的數目，在均勻選定a—"臓=30個角度{,，...,}，用多項式OJ擬合這些角度上的半徑{《,，&,}，使得所述每一層圖像的所述每一層輪廓線S表示為""p…，"c"V,、U"，…'A2"，用表示所有層圖像上/上得到的輪廓線，步驟(1.5).在所述所有層圖像/的輪廓線{5}中，按照高度均勻選定M層，以此表示所述所有層圖像/的外輪廓線，則相鄰兩層//,<//2之間，用下式表示任意一層輪廓線的坐標，用半徑i表示為<formula>formulaseeoriginaldocumentpage9</formula>，其中，A,A分別為在所述高度/Zp/^層角度a所對應的半徑，i^,/^分別是在高度巧,/^層所對應的12次多項式，步驟(1.6).把步驟(1.4)(1.5)得到的三維外輪廓線從極坐標[,凡//]轉換到三維笛卡爾坐標系[x,乂z]，得到了所述三維笛卡爾坐標系下的幾何模型G;步驟(2).把要進行活體熒光分子斷層成像的一個生物實驗用鼠放置在一個旋轉平臺上，激光器發出的激發光從一側照射所述活體生物實驗用鼠，激發光激發出熒光，被激發出的熒光在相對一側被CCD相機檢測，再把檢測結果輸入計算機中，分角度通過使所述旋轉平臺旋轉，采集不同旋轉角度被激發的熒光圖像，相當于實現了在不同角度時的激發光點光源激發；步驟(3).把歩驟(1)得到的幾何模型G導入一個有限元軟件包中，同時把所述各個激發光光源所對應的點加入所述幾何模型的節點上，再用有限元軟件包做網格劃分，離散成一個四面體網格Fem—ikfe^，該四面體網格的所有節點自動編號為!、1,…,W，iV是總的網格節點數；步驟(4).按以下步驟確定每個激發光點光源s所對應的在所述四面體網格i^m—M^/7上的有限元網格節點X說明書第4/14頁步驟(4.1).把垂直入射到所述活體生物實驗用鼠的激光束建模成入射光皮下一個自由傳輸自由程=1/(3/^(0)的位置r、.處的激發光點光源5(r-O，該5(r-O代表在位置rs,處的一個沖激函數，//,Ur、)是^處激發光波長的散射系數，歩驟(4.2).對于每個所述四面體網格Fem—7kfei上有限元網格節點/，若^,=。則該有限元網格節點/就是激發光點光源^所對應的有限元網格節點；步驟(5).確定所述每個激發光點光源s在所述Few—iWe^上有限元網格表面上的檢測點所對應的所述四面體網格Fem—上的有限元網格節點/:若所述表面節點/滿足以下關系^e{""g/e/<w""g7e—且a;e{/ze/g/z/1—/譜—，所述a"g7e—/，，"wg7e—/n.沐—/，,/e妙/1—/z妙為設定值，則所述節點/屬于所述激發光點光源^對應的網格節點；步驟(6).在得到步驟(5)所述各檢測點后，利用所述CCD相機檢測得到的熒光圖像，得到各檢測點發出的熒光圖像；步驟(7).針對所述給定位置^;處的激發光點光源5，按照以下歩驟計算對應所述激發光點光源s所述表面檢測點W1^2，…，t/W上檢測到的熒光強度同未知的熒光強度之間的線性關系，其中丄是檢測點總數步驟(7.1).按下式計算7VxiV維剛度矩陣《中的元素《[/,y]，該剛度矩陣&是生物組織在激發光波長在所述四面體網格i^w一7kfe^上的有限元節點上的吸收和散射系數的體現，該剛度矩陣《的行和列分別用,'和J'表示，所述的行/對應所述Fem—7kfe^上的所有節點按照編號順序排列而成的行，所述的列y'對應所述Fem一iWeA上的所有節點TV按照編號順序排列而成的列，下同，其中任意一元素《[/，/]表示為&["刀=J"。(A(>)▽W,0)Vy乂(/)+//"e0》A+由(r)^(r)A，其中，Q代表所述的幾何模型G內部對應的空間區域，3Q代表所述的幾何模型G的表面邊界，A(/)是位于H立置處激發光波長的漫射系數，所述馬(0=1/(3&)，；^(r)為位于r位置處的組織在所述激發光波長的散射系數，為已知值，/^(r)為位于r位置處的組織在所述激發光波長的吸收系數，為已知值，V,(。和^(r)是行編號為"列編號為y的有限元節點對應的Lagrange線性形函數，由所述有限元軟件包給出，VV,(r)=[3(r)/3y/,0)/,d(r)/3z]和Vy,(/)=[30)/5y,(/)/*，3y/,0)/&]分別是所述Lagrange線性形函數y,O)、^.(r)的梯度，《=0.1511，步驟(7.2).使用Matlab軟件包函數PCG按下式計算所述激發光點光源s在所述有限元網格Few—M;W節點上的強度分布向量^，(^是7Vxl列向《Ot'=S、，TVx1列向量5、的列編號為y的元素5,[y]用下式表示,fl,有限元節點_/位于所述;;.位置上，已知了在節點上的強度分布向量O,，得到任意位置^處的強度分步0￡(0=1]^[_/]步驟(7.3).按下式計算7VxN矩陣《的元素《[/,刀，該矩陣《是所述激發光強度0)￡(0對熒光的激發作用在所述四面體網格Few—MeW上的有限元節點上的體現，其中步驟(7.4).按下式計算WxiV維剛度矩陣i^中的元素i^[/,/1，該剛度矩陣i^表示生物組織在激發光波長在所述四面體網格Fe附—Me^上的有限元節點上的吸收和散射系數的體現，其中l['.，/]=J"。(化(。▽K0)'▽^0)+//。m(r》,(,(r》A+*J"3ny,(。&0)&'其中，DJ。是位于r位置處熒光波長的漫射系數，所述=1/(3/^0))，&(。為位于r位置處的組織在熒光波長的散射系數，為己知值，/^(r)為位于H立置處的組織在熒光波長的吸收系數，為已知值，步驟(7.5).假想當把一個熒光波長的點光源放在所述檢測點W處，t//eWl,"2,…,d"，使用Matlab軟件包函數PCG計算所產生的熒光場在所有所述四面體網格尸em—il^A上的有限元節點上的強度分布向量CE^，向量(D』是iVxl列向量，其中"'u」否則步驟(7.6).按照下式生成所述各檢測點^1，J2,…，^^處檢測到的熒光強度同未知的熒光探針分布之間的線性關系，用矩陣R表示步驟(8).針對空間順序分布的激發光點光源W,A^1,…,A^，以及相對應的所述檢測點，h=^，ivxi列向量a,的列編號為y的元素a,l/]用下式表示:按照下式構建總的線性關系矩陣『，下述。已知值.是各激發光點光源對應的被檢測熒光，為<formula>formulaseeoriginaldocumentpage12</formula>成Wx1列向量是未知的熒光探針在所述有限元網格Few—M^/z節點上的分布，Om,_是Mxl列向量，M所有激發光點光源對應的總的檢測點數目，『是Mx7V矩陣；步驟(9).按下述代數表達式，根據步驟(8)得到的總的線性關系矩陣『求解位置的熒光探針分布f/:<formula>formulaseeoriginaldocumentpage12</formula>分別是f/中第)個元素在第"—/f"+l和"—!'^次迭代時的值，其中<formula>formulaseeoriginaldocumentpage12</formula>t/[a'，是中第^個元素在第"次迭代時的值，是矩陣^第^行，第^列的值，『[/,，/1是矩陣W第^行，第y列的值，。,畫,用7表示，^,]是的r第A個元素，/1=0.1，C7迭代初始值設為^/=0。所述的朋g/e—/wv=-48°，""g/e—/z/g/z=48°，所述的/z".g/z/1—=-1.5c附,/7"'g似—/妙=1.5cw。為了驗證本重建方法的可靠性，實驗選擇了一套在實驗用鼠圖譜上得到的仿真數據。從重建結果可以看出，重建的熒光體具有良好的定位和定量效果。圖l熒光分子成像的機理l代表的是熒光濾光片。圖2熒光分子斷層成像的機理。圖3實驗用鼠圖譜的外表面輪廓。圖4非接觸熒光分子斷層成像系統結構圖l代表的是熒光濾光片。圖5實驗用鼠外表面輪廓劃分得到的有限元網格的部分顯示。圖6實驗用鼠圖譜模型上的仿真數據重建結果。圖7計算機程序的流程框圖。具體實施方式熒光分子成像的機理可以用圖1說明激發光光源發出的激發光激發組織內的熒光探針(標記的蛋白，腫瘤，熒光染料等)發射熒光，傳播到表面的熒光被CCD照相機檢測到。在圖l中，位于位置r、,的激發光點光源發出的激發光在實驗用鼠體內傳播，用實線箭頭表示。激發光遇到熒光探針后，會被熒光探針吸收，然后熒光探針被激發，發射出更高波長的熒光。熒光在實驗用鼠體內傳播，用虛線箭頭表示。到達激發光源對側實驗用鼠表面的熒光可以用CCD照相機檢測到。CCD照相機前有熒光濾光片，保證只有被激發的熒光才被CCD檢測到。由于生物組織具有強散射特性，因而激發光和被激發的熒光在實驗用鼠體內都不是直線傳播。這個成像的機理以及光在生物組織中的傳播可以用兩個耦合的稱為漫射方程的偏微分方程組來表示其中下標e代表激發光波長，下標m代表熒光波長。①e(。和OJr)分別代表激發光和被激發熒光的強度分布。r代表整個幾何模型內的空間位置坐標。在本文中，(...)表示函數操作，而[…]表示向量、圖像、或者矩陣操作。例如，(D。(0表示激發光光強場在r處的值，而o丄/]則代表激發光強場向量c^的第y個元素。光在組織中的傳播受到兩個因素的影響，稱為吸收和散射。吸收特性會吸收光的一部分能量，而散射特性會改變光的傳播方向，使得光在組織中的不能直線傳播。組織的吸收和散射特性在不同波長會不相同。A(。和Z)Jr)分別是位于r處的激發光和熒光波長的漫射系數，對應為D。(r)=l/(3Mre(0)和A>)=l/(3^>))，其中&(r)和&(0分別是組織在激發光和熒光波長的是散射系數。圪6("和/^一)分別是組織在激發光和熒光波長的吸收系數。在熒光分子斷層成像中，把垂直入射到實驗用鼠表面的激光束建模成入射點皮下一個傳輸自由程^=1/(3//〕處的位置r、.的激發光點光源50-O。-O代表在位置。處的一個沖激函數。t/(r)是位于r處的熒光探針的分布。熒光分子斷層成像是利用已知的在表面位置處檢測到的熒光強度，來重建未知的熒光探針的濃度分布C/(r)。其原理用圖2簡要說明分布在實驗用鼠體內的所有熒光探針對檢測點"l和d2處的檢測到的熒光強度都有貢獻，但是貢獻的權重會有不一樣。在圖中，用熒光探針1，2，3簡要代表實驗用鼠體內的所有熒光探針分布。以檢測點^1為例，熒光探針1，2，3對"l貢獻權重分別為f/lxwl，t/2xw2，t/3xw3。其中/71，f/2，f/3分別是熒光探針l，2，3的濃度，而wl，w2，w3是位于熒光探針l，2，3位置處的單位濃度的熒光探針對檢測點dl的貢獻權重，稱為單位貢獻權重。單位貢獻權重M，w2，w3可以通過熒光分子斷層重建方法求解得到。對不同的激發光源-熒光檢測點對，可以計算得到相應的單位貢獻權重。利用所有計算得到的單位貢獻權重以及檢測到的熒光強度，就可以重建出未知的熒光濃度f/1，f/2，t/3。這個單位貢獻權重在本文中被稱為檢測到的熒光強度同未知的熒光探針分布之間的線性關系。實際的重建方法計算有限元網格節點上的單位熒光濃度的熒光探針對檢測到的熒光強度的單位貢獻權重，然后可以用優化方法重建得到有限元網格節點上的熒光探針濃度分布。在熒光分子斷層成像中，多個角度的激發光點光源順序激發，同時檢測到相應的熒光投影數據，從多個角度投影數據就可以斷層重建出未知的熒光探針濃度分布。在熒光分子斷層成像系統中，用CT成像設備先得到實驗用鼠的解剖結構信息，然后利用醫學圖像分割的方法分割CT圖像，得到實驗用鼠的器官圖譜。這里提出的方法可以在得到的實驗用鼠圖譜的基礎上，利用測量到的多個角度的熒光投影數據，重建出未知的熒光探針濃度分布。方法的主要步驟包括步驟(l)，利用分割完畢的實驗用鼠圖譜，通過圖像處理得到外表面輪廓，并將外表面輪廓建模成幾何體，這個方法依次按以下步驟實現步驟(l.l).在已經得到的實驗用鼠圖譜中，圖譜中組織的灰度值都是大于等于1，而背景的灰度值是0。因而，選定一個閾值^"=0.5，就能夠分開圖譜中的組織和背景。對圖譜的每一層圖像/，利用閾值Ar進行二值化，大于閾值的取l，小于閾值的取0，得到二值圖5/:其中ww代表圖像的行，而"/代表圖像的列。步驟(1.2).用二值圖輪廓線提取的方法得到二值圖S/的外輪廓線S。方法如下對二值圖3/的每一個像素點A/[tow,co/]:(1)如果5/[raw,co/]為1，而且其相鄰的8個像素點至少有一個為0，那么這個像素點在外輪廓線S上；(2)如果6/[raw，co/]為0，而且其相鄰的8個像素點至少有一個為1，那么這個像素點在外輪廓線S上。這樣就得到了圖譜中一層圖像的外輪廓線。對圖譜中所有的層都執行二值化，然后進行輪廓線提取，就能得到所有層的外輪廓線，用{5}表示。步驟(1.3).依次把外輪廓線5*從像素坐標對應到笛卡爾坐標，給這些輪廓線選定一個中心[tow。，co/。]，每個像素點的幾何尺寸是Acm.那么以[row。,co/。]為原點，輪廓線上的任意一個點[tow,co/]都對應到幾何坐標[xj]:_y=(co/—co/0)*A然后把笛卡爾坐標[x，y]對應到極坐標[a,i]，其中《是極坐標系中的角度，A是極坐標系中的半徑。對每一條外輪廓線S上的所有點，每個點的半徑R用角度"按下式做12次多項式擬合凡二Za《，其中c為擬合次數，c=0,l，".12。a是擬合得到的多項式系數，用{&}表示。"=l,2,...Smw7，其中S—m/m是輪廓線S上的所有點的數目。在均勻選定a—mw^30個角度&,…,aa—用多項式{/^}擬合在這些角度上的半徑{^,'..,&■}。這樣，每一層圖像/所對應的輪廓線S表示為U",，…,"。—…A—步驟(1.4).每一層圖像都能夠得到一條輪廓線S，用{5}表示所有層圖像所得到的輪廓線。在{5}中，按照高度均勻選定M層，用這M層輪廓線來代表整個外輪廓。在相鄰兩層之間，可以通過插值的方法得到輪廓坐標。也就是說任意指定一個角度a和高度i/，都能夠相應計算得到其半徑i。設高度//處于高度為//,和//2的兩層(//1<//2)之間，那么其中/,,&分別是在Hp/7,層角度a所對應的半徑，/^和i^分別是在巧，A所對應的12次多項式。步驟(1.5).得到的外輪廓從三維極坐標[a,i，/Z]轉換到三維笛卡爾坐標[x，y,z]，這樣就得到了在笛卡爾坐標系下的幾何模型G，如圖3。步驟(2).在圖4所示的成像系統中，實驗用鼠被放置在一個旋轉平臺上。激光從一側照射實驗用鼠，激發實驗用鼠體內的熒光，被激發出的熒光在對側被CCD照相機檢測到。通過旋轉平臺，然后采集圖像，就能在不同角度實現激發，也就對應多個激發光點光源。把步驟(1)得到的幾何模型G導入有限元軟件包中，同時把各個激發光點光源所對應的點加到這個幾何模型G的節點上，然后用有限元軟件做網格劃分，離散成四面體網格i^附—MeA。圖5部分顯示了該四面體網格。該四面體網格的所有節點自動編號為/=1,一,^，iV是總的網格節點數。步驟(3).確定每個激發光點光源s以及其所有檢測點所對應的有限元節點。在離散的四面體網格&w一ikfei中，每個激發光點光源都對應一個有限元節點，通過以下方式確定激發光點光源對應的節點對每個有限元節點/，如果^、.=。那么/就是激發光點光源s所對應的有限元節點，其中^是激發光點光源s所對應的位置坐標，而r是有限元節點/所對應的位置坐標。同時，對每個激發光點光源5，在其對側對應的角度""g/e—/owtwg/e—Wg/以及高度/^'g似一/ow~/ze/g&—Wg/7范圍內檢測到熒光。在方法中，角度范圍選為-48'~48°，而高度范圍選為相對激發光點光源高度-1.51.5cm內。在方法中，要求檢測點放置在有限元網格的表面節點，因而，可以用以下方式確定檢測點對應的所有有限元節點對每一個表面節點/，如果ce{""g/e—/cw""g/e—/n.g/z}而且r,.e{/ze/g/^—/ow~—/z/g/}，另卩么節點/屬于這個光源對應的檢測點，否則不屬于。對特定激發光點光源，找到相應檢測點后，利用拍攝到的熒光圖象，得到在這些檢測點處的測量到的熒光強度。步驟(4).描述熒光分子斷層成像中光傳播的漫射方程組在有限元網格Fem—MeA上可以離散成一對線性方程組<formula>formulaseeoriginaldocumentpage16</formula>其中Wx1列向量cDt,和iVx1列向量。分別是是激發光和被激發熒光在網格節點上的強度分布；Wxl列向量"代表了熒光探針在網格節點上的分布；A^l列向量5、,代表了激發光點光源5(r-r)在網格節點的分布。在有限元方法中，己知了在離散的網格節點上的分布，可以通過有限元節點所對應的Lagmnge線性形函數來得到在任意位置r上的分布。以激發光強度為例，在任意位置r上的激發光強度(^(。，可以表示為0)=^0力']^(0，其中①』刀是列向量A的第7'個元素。①jy]也是激發光強度在編號為y的網格節點上的值。^,(r)是編號為的網格節點所對應的Lagrange線性形函數。&和《m在有限元方法中稱為剛度矩陣。矩陣^是組織在激發光波長的吸收和散射特性在網格節點上的體現。矩陣《是組織在熒光波長的吸收和散射特性在網格節點上的體現。矩陣《是激發光強度分布OJ。對熒光探針的激發作用在網格節點上的體現。這三個矩陣都是WxiV的稀疏矩陣。對于矩陣《，/^和《，它們的行和列分別用/和y表示，所述的行/對應Fem—T^fey/z上的所有節點7V按照編號順序排列而成的行，所述的列_/對應Few_ikfe^上的所有節點W按照編號順序排列而成的列。針對位于位置^;的激發光源s,其表面檢測節點WU2,…,^^上檢測到的熒光強度同未知的熒光探針分布之間的線性關系可以按照以下步驟得到步驟(4.1).計算剛度矩陣&。^的第/行，第y列元素&[/j]可以用下式求得，其中Q代表了整個幾何模型G內部所對應的空間區域；3Q則代表幾何模型G的表面邊界；下標e代表激發光波長；r是位置坐標，A(r)是位于r處的漫射系數，對應為化0)=1/(3&(r》；《和組織和空氣的折射系數有關，一般取g=0.1511;y,0)和^0)是行編號為/，列編號為j'的有限元節點對應的Lagrange線性形函數，由所述有限元軟件包給出。Vw,(r)=[3y,(r)/ax,3y,0)/^y,3y,(/)/&]和V(^(r)=[3<//,0)/&，3(^,0)//&]分別是所述Lagrange線性形函數^,(/)、^.(/)的梯度。稀疏矩陣《通過在有限元軟件包中設置好組織在激發光波長所對應的光學參數z^,00和mUO的接口函數，利用有限元軟件包求解得到。光學參數/^(r)和a。O0是這樣獲取的首先，獲取位置r的圖譜的灰度值，確定位置/"處所屬的組織器官；然后，查表得到該組織器官在激發光波長對應的光學參數。下表就是不同組織器官在激發光和熒光波長的光學參數值<table>tableseeoriginaldocumentpage18</column></row><table>其中，光學參數的單位是cnT1。歩驟(4.2).計算激發光點光源s所對應的激發光在有限元網格Few—ikfe^節點上的強度分布Ot,，《A=A。^的第_/個元素0)丄/]代表了激發光光強在編號為./的有限元節點上的值。列向量A的第./個元素5、[/l按下式表示，p,有限元節點y位于所述/;位置上力]^0，否則°利用雅克比預處理的共軛梯度法來迭代求解這個線性方程組，使用Matlab軟件包函數PCG來實現，輸入矩陣^和向量i,，輸出向量d^。步驟(4.3).計算矩陣F、。《的第z'行，第J'列元素《[/，)]表示為:其中ojo是激發光點光源^的光強場在空間位置r處的值，表示為o=^>。[刀1//，這個累加計算可以由有限元軟件包來完成。F、通過在有限元軟件包中設置好接口函數，利用有限元軟件包求解得到。這個接口函數是這樣的，輸入位置r，輸出r處的激發光強強度，。步驟(4.4).計算剛度矩陣《。矩陣尺第/行，第y列元素^[/，刀表示為，i["y]=J"q(ao)vw,(。.vw乂o)+(r)w,o》&+*j^y,o，乂。其中，下標m代表被激發光波長。這個矩陣的求解類似這個計算矩陣《，通過在有限元軟件包中設置好組織在被激發光波長的光學參數^m(r)和/^(r)接口函數，利用有限元軟件求解得到，光學參數Am(0和&(0是這樣獲取的首先，判斷位置r處所屬的組織器官；然后，査表得到該器官在熒光波長對應的光學參數。步驟(4.5).假設把一個熒光波長的點光源放在檢測節點^^/€{^1^2,...,^}處，計算產生的熒光場在有限元節點上的分布向量a^，其中iVxl列向量Sw的第J個元素5J/^'":n|。這個線性方程組通過雅克比預處理的共wJL0，否貝U軛梯度法來迭代求解，見步驟(4.2).步驟(4.6).生成檢測點{^/1^2,...,"}處檢測到的熒光同未知的熒光探針分布之間的線性關系矩陣『、.，可以看到，對于每個激發光源-熒光檢測對，例如卜J1，在檢測點dl檢測到的熒光同未知的熒光探針分布之間的線性關系通過一個向量-稀疏矩陣乘積。a尸,得到。步驟(5).針對多個順序激發的光源AJ=l...iV,以及其相應的檢測點，疊加生成總的線性關系矩陣『。按照歩驟(4)，分別建立每個激發光源A和其相應檢測點的線性關系矩陣『、.p然后把這些矩陣疊加可以得到總的線性關系矩陣，如下所示①■si其中，iVxl列向量t/是未知的熒光探針在有限元網格Few—7kfeW節點上的分布。Mxl列向量0),,。,是所有檢測點所對應的熒光強度值，其中M所有激發光點光源對應的總的檢測點數目。『是MxiV矩陣。步驟(6).根據線性關系矩陣重建未知的熒光探針分布，利用代數重建法(ART)來實現重建，把0^,，用向量r表示，求解線性方程組『17=「時，按照下式迭代<formula>formulaseeoriginaldocumentpage20</formula>其中"w、t/[y]",是;y中第乂個元素在第"—/^+i和"」^次迭代時的值。是t/中第^個元素在第"—/^次迭代時的值。嗎^,g是矩陣r第f,行，第^列的值。是矩陣ff第^行，第乂列的值。r[U是F第6個元素。義取值為義=0.1.迭代的初值取[/=0。示例的重建結果如圖6所示。圖6是重建的熒光探針分布在激發光點光源那個高度層上，坐標x,y各點畫出的等濃度線分布。權利要求1.一種基于線性關系的熒光分子斷層重建方法，其特征在于，所述方法是用一臺計算機依次按以下步驟實現的步驟(1).把用CT成像設備分割完畢的生物實驗用鼠的分層的圖譜輸入所述計算機中，該計算機通過下述圖像處理方法得到外表面輪廓，并把所表述表面輪廓建模成幾何體步驟(1.1).對所述圖譜的每一層圖像I，用一個相同的閾值thr＝0.5進行二值化大于所選閾值thr者取1，小于或等于該閾值thr者取0，把該圖譜中的組織和背景依次用所述的“1”和“0”分開，得到每一層圖像的二值圖BI其中，row代表該圖像I的行，col代表該圖像I的列，步驟(1.2).用二值圖輪廓線提取的方法得到所述的二值圖BI的外輪廓線S對該二值圖BI上的每一個像素點BI[row，col]執行以下操作若BI[row，col]為1，而且其相鄰的8個像素點至少有一個為0，則該像素點BI[row，col]在所述外輪廓線S上，若BI[row，col]為0，而且其相鄰的8個像素點至少有一個為1，則該像素點BI[row，col]在所述外輪廓線S上，步驟(1.3).把步驟(1.2)得到的所述外輪廓線S從像素坐標對應到笛卡爾坐標，則該外輪廓線S上的任意一個點[row，col]對應到所述笛卡爾坐標[x，y]上為，x＝(row-row0)*Δ，y＝(col-col0)*Δ，其中，點[row0，col0]為該輪廓線的中心，Δ是每個像素點的幾何尺寸，用cm表示，步驟(1.4).把步驟(1.3)得到的該笛卡爾坐標[x，y]對應到極坐標[α，R]上，α是角度，R是極坐標的半徑，對所述外輪廓線S上的所有點，把每個點上的所述半徑R用角度α按下式做多項式擬合c為擬合次數，c＝0，1，…12，pc是擬合得到的多項式系數，用{pc}表示，n＝1，2，…S_num，S_num是所述輪廓線S上的所有點的數目，在均勻選定α_num＝30個角度{α1，…，αα_num}，用多項式{pc}擬合這些角度上的半徑{R1，…，Rα_num}，使得所述每一層圖像的所述每一層輪廓線S表示為{{α1，…，αα_num}，{R1，…，Rα_num}，{p0，…，p12}}，用{S}表示所有層圖像上I上得到的輪廓線，步驟(1.5)在所述所有層圖像I的輪廓線{S}中，按照高度均勻選定M層，以此表示所述所有層圖像I的外輪廓線，則相鄰兩層H1<H2之間，用下式表示任意一層輪廓線的坐標，用半徑R表示為H1<H2，R＝R1*(H2-H)/(H2-H1)+R2(H-H1)/(H2-H1)，其中，R1，R2分別為在所述高度H1，H2層角度α所對應的半徑，PH1，RH2分別是在高度H1，H2層所對應的12次多項式，步驟(1.6).把步驟(1.4)～(1.5)得到的三維外輪廓線從極坐標[α，R，H]轉換到三維笛卡爾坐標系[x，y，z]，得到了所述三維笛卡爾坐標系下的幾何模型G；步驟(2).把要進行活體熒光分子斷層成像的一個生物實驗用鼠放置在一個旋轉平臺上，激光器發出的激發光從一側照射所述活體生物實驗用鼠，激發光激發出熒光，被激發出的熒光在相對一側被CCD相機檢測，再把檢測結果輸入計算機中，分角度通過使所述旋轉平臺旋轉，采集不同旋轉角度被激發的熒光圖像，相當于實現了在不同角度時的激發光點光源激發；步驟(3).把步驟(1)得到的幾何模型G導入一個有限元軟件包中，同時把所述各個激發光光源所對應的點加入所述幾何模型的節點上，再用有限元軟件包做網格劃分，離散成一個四面體網格Fem_Mesh，該四面體網格的所有節點自動編號為i＝1，…，N，N是總的網格節點數；步驟(4).按以下步驟確定每個激發光點光源s所對應的在所述四面體網格Fem_Mesh上的有限元網格節點步驟(4.1).把垂直入射到所述活體生物實驗用鼠的激光束建模成入射光皮下一個自由傳輸自由程的位置rs處的激發光點光源δ(r-rs)，該δ(r-rs)代表在位置rs處的一個沖激函數，是rs處激發光波長的散射系數，步驟(4.2).對于每個所述四面體網格Fem_Mesh上有限元網格節點i，若rs＝ri，則該有限元網格節點i就是激發光點光源s所對應的有限元網格節點；步驟(5).確定所述每個激發光點光源s在所述Fem_Mesh上有限元網格表面上的檢測點所對應的所述四面體網格Fem_Mesh上的有限元網格節點i若所述表面節點i滿足以下關系ri∈{angle_low～angle_high}且ri∈{height_low～height_high}，所述angle_low，angle_high，height_low，height_high為設定值，則所述節點i屬于所述激發光點光源s對應的網格節點；步驟(6).在得到步驟(5)所述各檢測點后，利用所述CCD相機檢測得到的熒光圖像，得到各檢測點發出的熒光圖像；步驟(7).針對所述給定位置rs處的激發光點光源s，按照以下步驟計算對應所述激發光點光源s所述表面檢測點{d1，d2，...，dL}上檢測到的熒光強度同未知的熒光強度之間的線性關系，其中L是檢測點總數步驟(7.1).按下式計算N×N維剛度矩陣Ke中的元素Ke[i，j]，該剛度矩陣Ke是生物組織在激發光波長在所述四面體網格Fem_Mesh上的有限元節點上的吸收和散射系數的體現，該剛度矩陣Ke的行和列分別用i和j表示，所述的行i對應所述Fem_Mesh上的所有節點N按照編號順序排列而成的行，所述的列j對應所述Fem_Mesh上的所有節點N按照編號順序排列而成的列，下同，其中任意一元素Ke[i，j]表示為其中，Ω代表所述的幾何模型G內部對應的空間區域，代表所述的幾何模型G的表面邊界，De(r)是位于r位置處激發光波長的漫射系數，所述為位于r位置處的組織在所述激發光波長的散射系數，為已知值，μae(r)為位于r位置處的組織在所述激發光波長的吸收系數，為已知值，ψi(r)和ψj(r)是行編號為i，列編號為j的有限元節點對應的Lagrange線性形函數，由所述有限元軟件包給出，和分別是所述Lagrange線性形函數ψi(r)、ψj(r)的梯度，q＝0.1511，步驟(7.2).使用Matlab軟件包函數PCG按下式計算所述激發光點光源s在所述有限元網格Fem_Mesh節點上的強度分布向量Φe，Φe是N×1列向量KeΦe＝Bs，N×1列向量Bs的列編號為j的元素Bs[j]用下式表示已知了在節點上的強度分布向量Φe，得到任意位置r處的強度分步步驟(7.3).按下式計算N×N矩陣Fs的元素Fs[i，j]，該矩陣Fs是所述激發光強度Φe(r)對熒光的激發作用在所述四面體網格Fem_Mesh上的有限元節點上的體現，其中Fs[i，j]＝∫ΩΦe(r)ψi(r)ψj(r)dr，步驟(7.4).按下式計算N×N維剛度矩陣Km中的元素Km[i，j]，該剛度矩陣Km表示生物組織在激發光波長在所述四面體網格Fem_Mesh上的有限元節點上的吸收和散射系數的體現，其中其中，Dm(r)是位于r位置處熒光波長的漫射系數，所述為位于r位置處的組織在熒光波長的散射系數，為已知值，μam(r)為位于r位置處的組織在熒光波長的吸收系數，為已知值，步驟(7.5).假想當把一個熒光波長的點光源放在所述檢測點dl處，dl∈{d1，d2，…，dL}，使用Matlab軟件包函數PCG計算所產生的熒光場在所有所述四面體網格Fem_Mesh上的有限元節點上的強度分布向量Φdl，向量Φdl是N×1列向量，其中KmΦdl＝Bdl，N×1列向量Bdl的列編號為j的元素Bdl[j]用下式表示步驟(7.6).按照下式生成所述各檢測點{d1，d2，…，dL}處檢測到的熒光強度同未知的熒光探針分布之間的線性關系，用矩陣Ws表示步驟(8).針對空間順序分布的激發光點光源sk，k＝1，…，Ns，以及相對應的所述檢測點，按照下式構建總的線性關系矩陣W，下述Φm，meas是各激發光點光源對應的被檢測熒光，為已知值N×1列向量U是未知的熒光探針在所述有限元網格Fem_Mesh節點上的分布，Φm，meas是M×1列向量，M所有激發光點光源對應的總的檢測點數目，W是M×N矩陣；步驟(9).按下述代數表達式，根據步驟(8)得到的總的線性關系矩陣W求解位置的熒光探針分布U其中U[j]tt_iter+1、U[j]tt_iter分別是U中第j個元素在第u_iter+1和u_iter次迭代時的值，U[t2]tt_iter是U中第t2個元素在第u_iter次迭代時的值，W[t1，t2]是矩陣W第t1行，第t2列的值，W[t1，j]是矩陣W第t1行，第j列的值，Φm，meas用V表示，V[t1]是的V第t1個元素，λ＝0.1，U迭代初始值設為U＝0。2."A'…'—…'W、"細〉，(/V…，Pi2"'用{S}表示所有層圖像上/上得到的輪廓線，步驟(1.5).在所述所有層圖像/的輪廓線{^中，按照高度均勻選定M層，以此表示所述所有層圖像/的外輪廓線，則相鄰兩層//,<//2之間，用下式表示任意一層輪廓線的坐標，用半徑i表示為W=《*(//2-77)/(//2—//,)+/2(//—/(//2—H》，其中，《,^2分別為在所述高度//1,//2層角度"所對應的半徑，/^,i^分別是在高度//,，//2層所對應的12次多項式，步驟(1.6).把步驟(1.4)(1.5)得到的三維外輪廓線從極坐標[a,i,/f]轉換到三維笛卡爾坐標系[;c,少,zl，得到了所述三維笛卡爾坐標系下的幾何模型G;步驟(2).把要進行活體熒光分子斷層成像的一個生物實驗用鼠放置在一個旋轉平臺上，激光器發出的激發光從一側照射所述活體生物實驗用鼠，激發光激發出熒光，被激發出的熒光在相對一側被CCD相機檢測，再把檢測結果輸入計算機中，分角度通過使所述旋轉平臺旋轉，采集不同旋轉角度被激發的熒光圖像，相當于實現了在不同角度時的激發光點光源激發；步驟(3).把步驟(1)得到的幾何模型G導入一個有限元軟件包中，同時把所述各個激發光光源所對應的點加入所述幾何模型的節點上，再用有限元軟件包做網格劃分，離散成一個四面體網格Few—ikfe^，該四面體網格的所有節點自動編號為z、l,…,iV，iV是總的網格節點數；步驟(4).按以下步驟確定每個激發光點光源s所對應的在所述四面體網格i^^^M^/z上的有限元網格節點步驟(4.1).把垂直入射到所述活體生物實驗用鼠的激光束建模成入射光皮下一個自由傳輸自由程=1/(3&(O)的位置r、.處的激發光點光源50-O，該-O代表在位置。處的一個沖激函數，/^(。是r、處激發光波長的散射系數，步驟(4.2).對于每個所述四面體網格&m—Ma/7上有限元網格節點/，若^=;，則該有限元網格節點/就是激發光點光源^所對應的有限元網格節點；步驟(5).確定所述每個激發光點光源s在所述Fem一Me^上有限元網格表面上的檢測點所對應的所述四面體網格Few—M^/z上的有限元網格節點/:3.則所述節點/屬于所述激發光點光源s對應的網格節點；歩驟(6).在得到歩驟(5)所述各檢測點后，利用所述CCD相機檢測得到的熒光圖像，得到各檢測點發出的熒光圖像；步驟(7).針對所述給定位置r、處的激發光點光源"按照以下步驟計算對應所述激發光點光源^所述表面檢測點{"1^2,...,^:}上檢測到的熒光強度同未知的熒光強度之間的線性關系，其中Z是檢測點總數步驟(7.1).按下式計算A^7V維剛度矩陣&中的元素A[/,/1，該剛度矩陣《是生物組織在激發光波長在所述四面體網格Few—MeA上的有限元節點上的吸收和散射系數的體現，該剛度矩陣《的行和列分別用/和_/表示，所述的行/對應所述&m—M^/z上的所有節點iV按照編號順序排列而成的行，所述的列/對應所述Fem—Afe^上的所有節點iV按照編號順序排列而成的列，下同，其中任意一元素《[/,y]表示為d["刀=J"q(AK(,)'▽0)+(0)^+由J"a^,(,V,，其中，Q代表所述的幾何模型G內部對應的空間區域，3Q代表所述的幾何模型G的表面邊界，A(。是位于r位置處激發光波長的漫射系數，所述AO)"/(3&(")，A。(。為位于H立置處的組織在所述激發光波長的散射系數，為已知值，/^(。為位于r位置處的組織在所述激發光波長的吸收系數，為已知值，y,(o和^.oo是行編號為/，列編號為y的有限元節點對應的Lagmnge線性形函數，由所述有限元軟件包給出，V^,.(r)=[3^,(。/&，3^/,(。/^,3^(。/&]禾口V=[5y/,(r)/分別是所述Lagrange線性形函數y/,0)、W,O)的梯度，^=0.1511,步驟(7.2).使用Matlab軟件包函數PCG按下式計算所述激發光點光源s在所述有限元網格Few—^fe^節點上的強度分布向量0。，d^是vVxl列向量已知了在節點上的強度分布向量^,,得到任意位置r處的強度分步Oe(。-Z氣[j']〖A=5、，Wxl列向量A的列編號為y的元素S、,[/l用下式表示nf.，fi,有限元節點y位于所述r,位置上，4.步驟(7.3).按下式計算WxN矩陣F、,的元素F、.[/,/1，該矩陣《是所述激發光強度^(r)對熒光的激發作用在所述四面體網格Few_Mm/2上的有限元節點上的體現，其中步驟(7.4).按下式計算A^W維剛度矩陣i^中的元素^J/,/1，該剛度矩陣/^表示生物組織在激發光波長在所述四面體網格/^w—M^/7上的有限元節點上的吸收和散射系數的體現，其中&['■，./]=J"。("m(。▽k0)▽+(Ok0》^+iJ"3。k0V,0)&'其中，DJO是位于r位置處熒光波長的漫射系數，所述^=1/(3^)，z^O)為位于r位置處的組織在熒光波長的散射系數，為已知值，為位于r位置處的組織在熒光波長的吸收系數，為已知值，步驟(7.5).假想當把一個熒光波長的點光源放在所述檢測點W處，^eW1^2,…^Z)，使用Matlab軟件包函數PCG計算所產生的熒光場在所有所述四面體網格i^w—7kfe^上的有限元節點上的強度分布向量0^，向量0^,是iVxl列向量，其中kA=A,，A^xl列向量^,的列編號為y'的元素A,[/l用下式表示，=1，7.=t//Lo，否則步驟(7.6).按照下式生成所述各檢測點Wl,d2,…,^^處檢測到的熒光強度同未知的熒光探針分布之間的線性關系，用矩陣R表示—①①F、,；步驟(8).針對空間順序分布的激發光點光源A,^:=1,...，乂，以及相對應的所述檢測點，按照下式構建總的線性關系矩陣『，下述o^,、是各激發光點光源對應的被檢測熒光，為已知值「l①.w、.iVxl列向量t/是未知的熒光探針在所述有限元網格Few7kfey/2節點上的分布，0>5.是Mxl列向量，M所有激發光點光源對應的總的檢測點數目，PF是MxiV矩陣；步驟(9).按下述代數表達式，根據步驟(8)得到的總的線性關系矩陣『求解位置的熒光探針分布"M廠w-iy^"2]"K]"-齡t/[乂]"--十1二"[乂]"-'旨+義￡^-『[qJ]射f/[刀"」"、"[刀"」""分別是C/中第/個元素在第"—+1和"—次迭代時的值，t/[g",是[/中第^個元素在第"—ter次迭代時的值，嗎^g是矩陣『第U亍，第^列的值，W[f,J]是矩陣W第f,行，第/列的值，氣，，用F表示，化]是的K第/,個元素，義=0.1，"迭代初始值設為^/=0。2.根據權利要求l所述的一種基于線性關系的熒光分子斷層重建方法，其特征在于，所述的""g/e—/ow=-48°，""g/e—/z/g//=48°，所述的/ze/g//"—/ow=-1.5cm,/ze/g/^—/z/g/=1.5c附。全文摘要一種基于線性關系的熒光分子斷層成像重建方法屬于在生物醫學中圖像重建的
技術領域：
。其特征在于，依次含有以下步驟一、實驗用鼠圖譜數據通過圖像處理和數值擬合得到外表面輪廓；二、這個外表面輪廓被導入有限元軟件包進行網格劃分；三、找出每個激發光源以及相應的熒光檢測點對應的有限元網格節點；四、針對一個特定的激發光源sk，建立表面檢測到的熒光強度同未知熒光探針分布之間的線性關系矩陣W_sk；五、針對所有激發光源，生成總的線性關系矩陣；六、得到這個線性關系后，代數重建方法被用來迭代逼近真實的熒光探針分布。本發明具有可處理非均勻光學參數以及快速建立表面檢測到的熒光強度同未知熒光探針分布之間的精確線性關系的特點。文檔編號G01N21/64GK101396262SQ20081022544公開日2009年4月1日申請日期2008年10月31日優先權日2008年10月31日發明者欣劉,李明澤,汪待發,王洪凱,凈白,燕金,陳延平,陳欣瀟申請人:清華大學