本發(fā)明涉及圖像處理技術(shù)領(lǐng)域，特別是涉及基于稀疏表示的改進(jìn)耦合字典學(xué)習(xí)的腦部ct/mr圖像融合方法。

背景技術(shù)：

在醫(yī)學(xué)領(lǐng)域，醫(yī)生需要對(duì)同時(shí)具有高空間和高光譜信息的單幅圖像進(jìn)行研究和分析，以便于對(duì)疾病進(jìn)行準(zhǔn)確診斷和治療。這種類型的信息僅從單模態(tài)圖像中無(wú)法獲取，例如，ct成像能夠捕捉人體的骨結(jié)構(gòu)，具有較高的分辨率，而mr成像能夠捕捉人體器官的軟組織如肌肉、軟骨、脂肪等細(xì)節(jié)信息。因此，將ct和mr圖像的互補(bǔ)信息相融合以獲取更全面豐富的圖像信息，可為臨床診斷和輔助治療提供有效幫助。

目前應(yīng)用于腦部醫(yī)學(xué)圖像融合領(lǐng)域比較經(jīng)典的方法是基于多尺度變換的方法：離散小波變換(dwt)、平穩(wěn)小波變換(swt)、雙樹復(fù)小波變換(dtcwt)、拉普拉斯金字塔(lp)、非下采樣contourlet變換(nsct)?；诙喑叨茸儞Q的方法能很好地提取圖像的顯著特征，但是對(duì)圖像誤配準(zhǔn)敏感，傳統(tǒng)的融合策略也使得融合結(jié)果無(wú)法保留源圖像的邊緣、紋理等細(xì)節(jié)信息。隨著壓縮感知的興起，基于稀疏表示的方法被廣泛用于圖像融合領(lǐng)域，并取得了極佳的融合效果。yang.b等使用冗余的dct字典稀疏表示源圖像，并使用“選擇最大”的規(guī)則融合稀疏系數(shù)。dct字典是一種由dct變換形成的隱式字典，易于快速實(shí)現(xiàn)，但是其表示能力有限。m.elad等提出k-svd算法用于從訓(xùn)練圖像中學(xué)習(xí)字典。與dct字典相比，學(xué)習(xí)的字典是一種自適應(yīng)于源圖像的顯式字典，具有較強(qiáng)的表示能力。學(xué)習(xí)的字典中將僅從自然圖像中采樣訓(xùn)練得到的字典稱為單字典，單字典可以表示任意一幅與訓(xùn)練樣本類別相同的自然圖像，但對(duì)于結(jié)構(gòu)復(fù)雜的腦部醫(yī)學(xué)圖像，使用單字典既要表示ct圖像又要表示mr圖像，難以得到精確的稀疏表示系數(shù)。b.ophir等提出小波域上的多尺度字典學(xué)習(xí)方法。即在小波域上對(duì)所有子帶分別使用k-svd算法訓(xùn)練得到所有子帶對(duì)應(yīng)的子字典。多尺度字典有效地將解析字典和學(xué)習(xí)的字典的優(yōu)勢(shì)相結(jié)合，能夠捕捉圖像在不同尺度和不同方向上包含的不同特征。但是所有子帶的子字典也是單字典，使用子字典對(duì)所有子帶進(jìn)行稀疏表示仍然難以得到精確的稀疏表示系數(shù)，并且分離的字典學(xué)習(xí)時(shí)間效率較低。yu.n等提出基于聯(lián)合稀疏表示的圖像融合方法兼具去噪功能。這種方法是對(duì)待融合的源圖像本身進(jìn)行字典學(xué)習(xí)，根據(jù)jsm-1模型提取待融合圖像的公共特征和各自特征，再分別組合并重構(gòu)得到融合圖像。這種方法由于是對(duì)待融合源圖像本身訓(xùn)練字典所以適用于腦部醫(yī)學(xué)圖像，可以得到精確的稀疏表示系數(shù)。但是對(duì)于每對(duì)待融合的源圖像都需要訓(xùn)練字典，時(shí)間效率低，缺乏靈活性。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明實(shí)施例提供了基于稀疏表示的改進(jìn)耦合字典學(xué)習(xí)的腦部ct/mr圖像融合方法，可以解決現(xiàn)有技術(shù)中存在的問題。

一種基于稀疏表示的改進(jìn)耦合字典學(xué)習(xí)的腦部ct/mr圖像融合方法，該方法包括：

預(yù)處理階段：對(duì)于已經(jīng)配準(zhǔn)的腦部ct/mr源圖像ic，ir∈rmn，rmn表示具有m行n列的向量空間，使用步長(zhǎng)為1的滑動(dòng)窗把源圖像ic，ir分別分割為大小的圖像塊，對(duì)于每幅ct源圖像ic和mr源圖像ir，都有個(gè)圖像塊，然后將這些圖像塊編纂成m維列向量，將ct源圖像ic中的第j個(gè)圖像塊記為mr源圖像ir中的第j個(gè)圖像塊記為減去各自的平均值：

其中，和分別表示和中所有元素的均值，1表示一個(gè)全1的m維列向量；

融合階段：使用coefromp算法求解的稀疏系數(shù)，公式表示如下：

其中，||α||0表示稀疏系數(shù)α中非零元素的個(gè)數(shù)，ε表示允許偏差的精度，d^f表示字典d^c和d^r融合后得到的融合字典；

將稀疏系數(shù)的l2范數(shù)作為源圖像的活躍度測(cè)量，則稀疏系數(shù)和通過以下的融合規(guī)則融合：

均值和使用“加權(quán)平均”規(guī)則融合：

其中，則和的融合結(jié)果為：

重建階段：對(duì)所有的圖像塊都執(zhí)行預(yù)處理階段和融合階段以得到所有圖像塊的融合結(jié)果，對(duì)于每個(gè)塊向量通過反滑窗的過程重塑成的圖像塊并放回到對(duì)應(yīng)的像素位置，再對(duì)重復(fù)像素取平均得到最終的融合圖像if。

優(yōu)選地，在融合階段中，所述融合字典通過以下方法計(jì)算獲得：

使用高質(zhì)量ct和mr圖像作為訓(xùn)練集，從訓(xùn)練集中采樣得到向量對(duì){x^c,x^r}，定義為n個(gè)采樣的ct圖像向量組成的矩陣，為對(duì)應(yīng)的n個(gè)采樣的mr圖像向量組成的矩陣，其中r^d×n表示具有d行n列的向量空間；

在字典學(xué)習(xí)代價(jià)函數(shù)基礎(chǔ)上加入支撐完整的先驗(yàn)信息，交替地更新d^c，d^r和a，對(duì)應(yīng)的訓(xùn)練優(yōu)化問題如下：

其中，a是x^c和x^r的聯(lián)合稀疏系數(shù)矩陣，τ是聯(lián)合稀疏系數(shù)矩陣a的稀疏度，⊙表示點(diǎn)乘，掩膜矩陣m由元素0和1組成，定義為m＝{|a|＝0}，等價(jià)于如果a(i,j)＝0則m(i,j)＝1，否則為0，引入輔助變量：

則式(1)可以等價(jià)的轉(zhuǎn)化為：

式(3)的求解過程分為稀疏編碼和字典更新兩個(gè)步驟：

首先，在稀疏編碼階段，隨機(jī)矩陣初始化字典和通過求解式(4)來(lái)實(shí)現(xiàn)對(duì)聯(lián)合稀疏系數(shù)矩陣a的更新：

分別對(duì)聯(lián)合稀疏系數(shù)矩陣a中每一列的非零元素進(jìn)行處理，而保持零元素完備，則式(4)可以轉(zhuǎn)換為下式：

式中，是對(duì)應(yīng)a的非零支集的子矩陣，αi是a第i列的非零部分，式(5)由系數(shù)重用正交匹配追蹤算法coefromp求解，得到更新的聯(lián)合稀疏系數(shù)矩陣a；

其次，在字典更新階段，式(3)的優(yōu)化問題轉(zhuǎn)化為：

則式(6)的補(bǔ)償項(xiàng)寫為：

式中，表示字典中待更新的第k列，表示聯(lián)合稀疏系數(shù)矩陣a的第k行，表示掩膜矩陣m的第j行，用來(lái)保證中的零元素在正確的位置，掩膜矩陣是將行向量復(fù)制d次得到尺寸為d×n的秩為1的矩陣，利用掩膜矩陣可以有效地去除中那些未用到第k個(gè)原子所對(duì)應(yīng)樣本的列，對(duì)誤差矩陣ek進(jìn)行奇異值分解(svd)得到ek＝uδv^t，使用矩陣u的第一列更新字典中的原子同時(shí)將稀疏系數(shù)矩陣a中的更新為矩陣v的第一列與δ(1,1)的乘積；

最后，循環(huán)執(zhí)行稀疏編碼和字典更新這兩個(gè)階段，直至達(dá)到預(yù)設(shè)的迭代次數(shù)為止，輸出一對(duì)耦合的d^c和d^r字典。

優(yōu)選地，使用以下方法對(duì)字典d^c和d^r進(jìn)行融合：

l^c(n)和l^r(n)，n＝1,2,…,n分別代表ct字典和mr字典的第n個(gè)原子的特征指標(biāo)，融合公式表示如下：

此處設(shè)λ＝0.25。

本發(fā)明實(shí)施例提供的基于稀疏表示的改進(jìn)耦合字典學(xué)習(xí)的腦部ct/mr圖像融合方法，能分別對(duì)正常腦部、腦萎縮和腦腫瘤三組腦部醫(yī)學(xué)圖像進(jìn)行融合，多次實(shí)驗(yàn)結(jié)果表明本發(fā)明提出的icdl方法與基于多尺度變換的方法、傳統(tǒng)稀疏表示的方法、基于k-svd字典學(xué)習(xí)的方法以及多尺度字典學(xué)習(xí)的方法相比，不僅提高了腦部醫(yī)學(xué)圖像融合的質(zhì)量，而且有效降低了字典訓(xùn)練的時(shí)間，能為臨床醫(yī)療診斷提供有效幫助。

附圖說明

為了更清楚地說明本發(fā)明實(shí)施例或現(xiàn)有技術(shù)中的技術(shù)方案，下面將對(duì)實(shí)施例或現(xiàn)有技術(shù)描述中所需要使用的附圖作簡(jiǎn)單地介紹，顯而易見地，下面描述中的附圖僅僅是本發(fā)明的一些實(shí)施例，對(duì)于本領(lǐng)域普通技術(shù)人員來(lái)講，在不付出創(chuàng)造性勞動(dòng)的前提下，還可以根據(jù)這些附圖獲得其他的附圖。

圖1為本發(fā)明實(shí)施例提供的基于稀疏表示的改進(jìn)耦合字典學(xué)習(xí)的腦部ct/mr圖像融合方法流程圖；

圖2為作為訓(xùn)練集的高質(zhì)量的ct和mr圖像；

圖3為正常腦部的ct/mr融合結(jié)果，其中a為ct圖像，b為mr圖像，c為dwt(離散小波變換)圖像，d為swt(平穩(wěn)小波變換)圖像，e為nsct(非下采樣contourlet變換)圖像，f為srm(傳統(tǒng)稀疏表示方法)圖像，g為srk(基于k-svd字典學(xué)習(xí)方法)圖像，h為mdl(基于多尺度字典學(xué)習(xí)方法)圖像，i為本發(fā)明使用的icdl(improvedcoupleddictionarylearning)圖像；

圖4為腦萎縮的ct/mr融合結(jié)果；

圖5為腦腫瘤的ct/mr融合結(jié)果。

具體實(shí)施方式

下面將結(jié)合本發(fā)明實(shí)施例中的附圖，對(duì)本發(fā)明實(shí)施例中的技術(shù)方案進(jìn)行清楚、完整地描述，顯然，所描述的實(shí)施例僅僅是本發(fā)明一部分實(shí)施例，而不是全部的實(shí)施例?；诒景l(fā)明中的實(shí)施例，本領(lǐng)域普通技術(shù)人員在沒有做出創(chuàng)造性勞動(dòng)前提下所獲得的所有其他實(shí)施例，都屬于本發(fā)明保護(hù)的范圍。

參照?qǐng)D1，本發(fā)明實(shí)施例中提供的基于稀疏表示的改進(jìn)耦合字典學(xué)習(xí)的腦部ct/mr圖像融合方法，該方法包括以下步驟：

步驟100，預(yù)處理階段：對(duì)于已經(jīng)配準(zhǔn)的腦部ct/mr源圖像ic，ir∈rmn，rmn表示具有m行n列的向量空間，使用步長(zhǎng)為1的滑動(dòng)窗把源圖像ic，ir分別分割為大小的圖像塊，對(duì)于每幅ct源圖像ic和mr源圖像ir，都有個(gè)圖像塊，然后將這些圖像塊編纂成m維列向量，將ct源圖像ic中的第j個(gè)圖像塊記為mr源圖像ir中的第j個(gè)圖像塊記為減去各自的平均值：

其中，和分別表示和中所有元素的均值，1表示一個(gè)全1的m維列向量；

步驟200，融合階段：使用coefromp算法求解的稀疏系數(shù)，公式表示如下：

其中，||α||0表示稀疏系數(shù)α中非零元素的個(gè)數(shù)，ε表示允許偏差的精度，d^f表示字典d^c和d^r融合后得到的融合字典，其具體計(jì)算方法如下：

使用圖2所示的高質(zhì)量ct和mr圖像作為訓(xùn)練集，從訓(xùn)練集中采樣得到向量對(duì){x^c,x^r}，定義為n個(gè)采樣的ct圖像向量組成的矩陣，為對(duì)應(yīng)的n個(gè)采樣的mr圖像向量組成的矩陣，其中r^d×n表示具有d行n列的向量空間；

本發(fā)明的耦合字典訓(xùn)練使用改進(jìn)的k-svd算法，該算法在傳統(tǒng)的字典學(xué)習(xí)代價(jià)函數(shù)基礎(chǔ)上加入支撐完整的先驗(yàn)信息，交替地更新d^c，d^r和a，對(duì)應(yīng)的訓(xùn)練優(yōu)化問題如下：

其中，a是x^c和x^r的聯(lián)合稀疏系數(shù)矩陣，τ是聯(lián)合稀疏系數(shù)矩陣a的稀疏度，⊙表示點(diǎn)乘，掩膜矩陣m由元素0和1組成，定義為m＝{|a|＝0}，等價(jià)于如果a(i,j)＝0則m(i,j)＝1，否則為0。因此a⊙m＝0能使a中所有零項(xiàng)保持完備。引入輔助變量：

則式(3)可以等價(jià)的轉(zhuǎn)化為：

式(5)的求解過程分為稀疏編碼和字典更新兩個(gè)步驟。

首先，在稀疏編碼階段，隨機(jī)矩陣初始化字典和通過求解式(6)來(lái)實(shí)現(xiàn)對(duì)聯(lián)合稀疏系數(shù)矩陣a的更新：

分別對(duì)聯(lián)合稀疏系數(shù)矩陣a中每一列的非零元素進(jìn)行處理，而保持零元素完備，則式(6)可以轉(zhuǎn)換為下式：

式中，是對(duì)應(yīng)a的非零支集的子矩陣，αi是a第i列的非零部分。式(7)由系數(shù)重用正交匹配追蹤算法coefromp求解，由此可得到更新的聯(lián)合稀疏系數(shù)矩陣a。

其次，在字典更新階段，式(5)的優(yōu)化問題可以轉(zhuǎn)化為：

則式(8)的補(bǔ)償項(xiàng)可以寫為：

式中，表示字典中待更新的第k列，表示聯(lián)合稀疏系數(shù)矩陣a的第k行，表示掩膜矩陣m的第j行，用來(lái)保證中的零元素在正確的位置。掩膜矩陣是將行向量復(fù)制d次得到尺寸為d×n的秩為1的矩陣，利用掩膜矩陣可以有效地去除中那些未用到第k個(gè)原子所對(duì)應(yīng)樣本的列。對(duì)誤差矩陣ek進(jìn)行奇異值分解(svd)得到ek＝uδv^t，使用矩陣u的第一列更新字典中的原子同時(shí)將稀疏系數(shù)矩陣a中的更新為矩陣v的第一列與δ(1,1)的乘積。

最后，循環(huán)執(zhí)行稀疏編碼和字典更新這兩個(gè)階段，直至達(dá)到預(yù)設(shè)的迭代次數(shù)為止，輸出一對(duì)耦合的d^c和d^r字典。然后使用以下方法對(duì)字典d^c和d^r進(jìn)行融合：

l^c(n)和l^r(n)，n＝1,2,…,n分別代表ct字典和mr字典的第n個(gè)原子的特征指標(biāo)，由于腦部ct和mr圖像是對(duì)應(yīng)于人體同一部位由不同成像設(shè)備獲取的圖像，因此兩者之間一定存在著公共特征和各自特征。本發(fā)明提出將特征指標(biāo)相差較大的原子看作各自特征，使用“選擇最大”規(guī)則融合。特征指標(biāo)相差較小的原子看作公共特征，使用“平均”的規(guī)則融合，公式表示如下：

此處設(shè)λ＝0.25，根據(jù)醫(yī)學(xué)圖像的物理特性，使用信息熵作為特征指標(biāo)。這種方法將稀疏域和空間域的方法結(jié)合起來(lái)，考慮了醫(yī)學(xué)圖像的物理特性計(jì)算字典原子的特征指標(biāo)，與稀疏域的方法相比，具有更加明確的物理意義。

由于字典更新階段同時(shí)更新字典和稀疏表示系數(shù)的非零元素，使得字典的表示誤差更小且字典的收斂速度更快。在稀疏編碼階段，考慮到每次迭代時(shí)都忽略前一次迭代的表示，coefromp算法提出利用上次迭代的稀疏表示殘差信息進(jìn)行系數(shù)更新，從而更快地得到所要求問題的解。

計(jì)算得到融合字典d^f后，將稀疏系數(shù)的l2范數(shù)作為源圖像的活躍度測(cè)量，則稀疏系數(shù)和通過以下的融合規(guī)則融合：

均值和使用“加權(quán)平均”規(guī)則融合：

其中，則和的融合結(jié)果為：

步驟300，重建階段：對(duì)所有的圖像塊都執(zhí)行上述兩個(gè)步驟以得到所有圖像塊的融合結(jié)果。對(duì)于每個(gè)塊向量通過反滑窗的過程重塑成的圖像塊并放回到對(duì)應(yīng)的像素位置，再對(duì)重復(fù)像素取平均得到最終的融合圖像if。

為驗(yàn)證本發(fā)明方法的有效性，選取三組已經(jīng)配準(zhǔn)的腦部ct/mr圖像進(jìn)行融合，分別為正常腦部ct/mr如圖3中a和b所示，腦萎縮ct/mr如圖4中a和b所示，腦腫瘤ct/mr如圖5中a和b所示，圖片大小均為256×256。選取的對(duì)比算法有：離散小波變換(dwt)、平穩(wěn)小波變換(swt)、非下采樣contourlet變換(nsct)、傳統(tǒng)稀疏表示的方法(srm)、基于k-svd字典學(xué)習(xí)的方法(srk)、基于多尺度字典學(xué)習(xí)的方法(mdl)，融合結(jié)果分別見圖3中c、d、e、f、g、h，圖4中c、d、e、f、g、h和圖5中c、d、e、f、g、h。

基于多尺度變換的方法中，對(duì)于dwt和swt方法，分解水平都設(shè)為3,小波基分別設(shè)為“db6”和“bior1.1”。nsct方法使用“9-7”金字塔濾波器和“c-d”方向?yàn)V波器，分解水平設(shè)為{2²,2²,2³,2⁴}?；谙∈璞硎镜姆椒ㄖ谢瑒?dòng)步長(zhǎng)為1，圖像塊大小均為8×8，字典大小均為64×256，誤差ε＝0.01，稀疏度τ＝6，icdl方法使用改進(jìn)的k-svd算法，執(zhí)行6個(gè)多重字典更新周期(duc)和30次迭代。

由圖3-5可以看出，dwt方法的融合圖像邊緣紋理模糊，圖像信息失真且存在塊效應(yīng)；與dwt方法相比，swt和nsct方法的融合質(zhì)量相對(duì)較好，圖像的亮度、對(duì)比度、清晰度有了很大的提升，但仍存在邊緣亮度失真，軟組織和病灶區(qū)域存在偽影的問題；srm和srk方法較基于多尺度變換的方法圖像的骨組織和軟組織更加清晰，偽影也有所減少，能很好地識(shí)別病灶區(qū)域；mdl方法與srm和srk方法相比能保留更多的細(xì)節(jié)信息，圖像質(zhì)量取得進(jìn)一步的改善，但仍有部分偽影存在；本發(fā)明提出的icdl方法在圖像的亮度、對(duì)比度、清晰度和細(xì)節(jié)的保持度上都優(yōu)于其他方法，融合圖像沒有偽影，骨組織、軟組織和病灶區(qū)域顯示清晰，有助于醫(yī)生診斷。

本領(lǐng)域內(nèi)的技術(shù)人員應(yīng)明白，本發(fā)明的實(shí)施例可提供為方法、系統(tǒng)、或計(jì)算機(jī)程序產(chǎn)品。因此，本發(fā)明可采用完全硬件實(shí)施例、完全軟件實(shí)施例、或結(jié)合軟件和硬件方面的實(shí)施例的形式。而且，本發(fā)明可采用在一個(gè)或多個(gè)其中包含有計(jì)算機(jī)可用程序代碼的計(jì)算機(jī)可用存儲(chǔ)介質(zhì)(包括但不限于磁盤存儲(chǔ)器、cd-rom、光學(xué)存儲(chǔ)器等)上實(shí)施的計(jì)算機(jī)程序產(chǎn)品的形式。

本發(fā)明是參照根據(jù)本發(fā)明實(shí)施例的方法、設(shè)備(系統(tǒng))、和計(jì)算機(jī)程序產(chǎn)品的流程圖和/或方框圖來(lái)描述的。應(yīng)理解可由計(jì)算機(jī)程序指令實(shí)現(xiàn)流程圖和/或方框圖中的每一流程和/或方框、以及流程圖和/或方框圖中的流程和/或方框的結(jié)合?？商峁┻@些計(jì)算機(jī)程序指令到通用計(jì)算機(jī)、專用計(jì)算機(jī)、嵌入式處理機(jī)或其他可編程數(shù)據(jù)處理設(shè)備的處理器以產(chǎn)生一個(gè)機(jī)器，使得通過計(jì)算機(jī)或其他可編程數(shù)據(jù)處理設(shè)備的處理器執(zhí)行的指令產(chǎn)生用于實(shí)現(xiàn)在流程圖一個(gè)流程或多個(gè)流程和/或方框圖一個(gè)方框或多個(gè)方框中指定的功能的裝置。

這些計(jì)算機(jī)程序指令也可存儲(chǔ)在能引導(dǎo)計(jì)算機(jī)或其他可編程數(shù)據(jù)處理設(shè)備以特定方式工作的計(jì)算機(jī)可讀存儲(chǔ)器中，使得存儲(chǔ)在該計(jì)算機(jī)可讀存儲(chǔ)器中的指令產(chǎn)生包括指令裝置的制造品，該指令裝置實(shí)現(xiàn)在流程圖一個(gè)流程或多個(gè)流程和/或方框圖一個(gè)方框或多個(gè)方框中指定的功能。

這些計(jì)算機(jī)程序指令也可裝載到計(jì)算機(jī)或其他可編程數(shù)據(jù)處理設(shè)備上，使得在計(jì)算機(jī)或其他可編程設(shè)備上執(zhí)行一系列操作步驟以產(chǎn)生計(jì)算機(jī)實(shí)現(xiàn)的處理，從而在計(jì)算機(jī)或其他可編程設(shè)備上執(zhí)行的指令提供用于實(shí)現(xiàn)在流程圖一個(gè)流程或多個(gè)流程和/或方框圖一個(gè)方框或多個(gè)方框中指定的功能的步驟。

盡管已描述了本發(fā)明的優(yōu)選實(shí)施例，但本領(lǐng)域內(nèi)的技術(shù)人員一旦得知了基本創(chuàng)造性概念，則可對(duì)這些實(shí)施例作出另外的變更和修改。所以，所附權(quán)利要求意欲解釋為包括優(yōu)選實(shí)施例以及落入本發(fā)明范圍的所有變更和修改。

顯然，本領(lǐng)域的技術(shù)人員可以對(duì)本發(fā)明進(jìn)行各種改動(dòng)和變型而不脫離本發(fā)明的精神和范圍。這樣，倘若本發(fā)明的這些修改和變型屬于本發(fā)明權(quán)利要求及其等同技術(shù)的范圍之內(nèi)，則本發(fā)明也意圖包含這些改動(dòng)和變型在內(nèi)。