寧 騫 徐 桓 劉 洋 高 鵬 劉本源 劉銳崗 尤富生 濮黃生*

基于雙波長(zhǎng)減影算法的X射線發(fā)光成像研究*

寧 騫①徐 桓②劉 洋①高 鵬①劉本源①劉銳崗①尤富生①濮黃生①*

目的：探討雙波長(zhǎng)減影算法對(duì)X射線發(fā)光成像(XLI)中存在的環(huán)境噪聲、X射線噪聲以及非目標(biāo)熒光信號(hào)等背景噪聲的去除效果。方法：開(kāi)發(fā)基于雙波長(zhǎng)減影算法的XLI圖像處理方法，并通過(guò)數(shù)值仿真以及仿體實(shí)驗(yàn)來(lái)對(duì)所提方法進(jìn)行驗(yàn)證。使用硬閾值算法和改進(jìn)型中值濾波算法對(duì)原始圖像進(jìn)行預(yù)處理，然后使用雙波長(zhǎng)減影算法處理，并對(duì)處理結(jié)果進(jìn)行評(píng)估。結(jié)果：在數(shù)值仿真中XLI圖像的目標(biāo)背景比(TBR)提高約38倍，在仿體實(shí)驗(yàn)中TBR提高約4倍。數(shù)值仿真以及仿體實(shí)驗(yàn)表明，所提方法處理圖像之后可以提升XLI圖像的TBR，改善圖像質(zhì)量。結(jié)論：基于雙波長(zhǎng)減影算法的圖像處理方法可以有效去除XLI圖像中的背景噪聲干擾，提取目標(biāo)熒光信號(hào)。

X射線發(fā)光成像；光學(xué)分子成像；背景噪聲；雙波長(zhǎng)減影；目標(biāo)背景比

X射線發(fā)光成像(X-ray luminescence imaging，XLI)是一種新興的應(yīng)用X射線激發(fā)稀土納米探針發(fā)光的光學(xué)分子成像技術(shù)，而稀土納米探針是一種摻雜稀土元素的磷納米顆粒[1-3]。當(dāng)使用X射線照射時(shí)，磷納米顆粒會(huì)發(fā)出可以穿透生物組織并被光電探測(cè)設(shè)備探測(cè)到的光信號(hào)，從而可對(duì)成像目標(biāo)進(jìn)行成像。XLI相較于現(xiàn)有的分子影像技術(shù)具有拓展熒光分子探針的合成種類(lèi)和應(yīng)用范圍、結(jié)構(gòu)與功能的雙模態(tài)成像以及自體熒光對(duì)成像質(zhì)量的影響小的優(yōu)勢(shì)[4-7]。但是，XLI在成像過(guò)程中通常會(huì)受到背景噪聲影響，其背景噪聲包括環(huán)境噪聲、X射線噪聲以及非目標(biāo)熒光信號(hào)等。環(huán)境噪聲主要為電荷耦合器件(charged coupled device，CCD)相機(jī)的本底噪聲和環(huán)境光造成的噪聲等；X射線噪聲指使用X射線激發(fā)稀土納米探針時(shí)泄露的X射線造成的高強(qiáng)度噪聲；非目標(biāo)熒光信號(hào)包括其他探針熒光、二次自體熒光等熒光信號(hào)。背景噪聲會(huì)降低圖像的對(duì)比度，使目標(biāo)熒光信號(hào)不易區(qū)分，甚至完全掩蓋目標(biāo)信號(hào)，從而降低成像結(jié)果的可識(shí)別度和可信度。

圖像減影算法是一種利用不同圖像之間差異來(lái)進(jìn)行感興趣目標(biāo)提取的算法，可以簡(jiǎn)單、有效地從背景噪聲信號(hào)中提取出目標(biāo)成像信號(hào)[8-9]。稀土納米探針的發(fā)光譜主要由摻雜的稀土元素決定，不同的稀土元素具有不同的發(fā)光特性，如摻雜銪(Eu)的磷納米顆粒發(fā)射波峰主要在紅光波段，而摻雜鋱(Tb)的磷納米顆粒發(fā)射波峰主要在綠光波段[10-11]。本研究基于不同磷納米顆粒的發(fā)射光譜差異，提出了一種雙波長(zhǎng)圖像減影算法來(lái)進(jìn)行背景噪聲濾除，進(jìn)而提取出目標(biāo)熒光信號(hào)。

1 雙波長(zhǎng)圖像減影算法

該部分為圖像處理所涉及的算法及其原理介紹，主要包括硬閾值算法、改進(jìn)型中值濾波算法以及雙波長(zhǎng)減影算法3個(gè)部分。算法的開(kāi)發(fā)及編程均在Matlab2011語(yǔ)言環(huán)境下進(jìn)行。

1.1 硬閾值算法

針對(duì)環(huán)境噪聲強(qiáng)度較低的特點(diǎn)，本研究擬采取硬閾值算法去除環(huán)境噪聲。首先確定背景區(qū)域，然后求取背景區(qū)域內(nèi)所有像素點(diǎn)數(shù)值的平均值m作為當(dāng)前實(shí)驗(yàn)條件下環(huán)境噪聲閾值；最后對(duì)整幅圖像進(jìn)行處理，假設(shè)圖像中像素點(diǎn)數(shù)值為k，如果k≤m，則賦值k=0；否則，k=k-m。

1.2 改進(jìn)型中值濾波算法

針對(duì)X射線噪聲高強(qiáng)度及稀疏特性，本研究采用中值濾波算法去除X射線噪聲，但是如果對(duì)整幅圖像進(jìn)行濾波處理，容易將圖像中的邊緣細(xì)節(jié)誤認(rèn)為噪聲點(diǎn)，可能會(huì)造成圖像中非噪聲區(qū)域失真。為了實(shí)現(xiàn)盡可能保持圖像中細(xì)節(jié)信息的同時(shí)去除高強(qiáng)度噪聲，本研究提出一種基于閾值篩選的改進(jìn)型中值濾波算法。首先采用閾值法對(duì)X射線噪聲點(diǎn)進(jìn)行篩選，然后用中值濾波算法對(duì)篩選出的噪聲點(diǎn)進(jìn)行賦值。

(1)X射線噪聲點(diǎn)篩選(閾值法)。X射線噪聲點(diǎn)篩選主要難點(diǎn)在于閾值的選擇，依據(jù)先驗(yàn)知識(shí)可知，X射線噪聲點(diǎn)個(gè)數(shù)不超過(guò)整幅圖像像素點(diǎn)個(gè)數(shù)的10%，且X射線噪聲點(diǎn)強(qiáng)度M(noi)為感興趣區(qū)(region of interest，ROI)內(nèi)非噪聲信號(hào)強(qiáng)度平均值M(sig)的k(k∈[5，10])倍，其計(jì)算為公式1：

可求得X射線噪聲點(diǎn)強(qiáng)度M(noi)與ROI內(nèi)信號(hào)強(qiáng)度平均值M(ROI)之比k(nr)為公式2：

因此，本次實(shí)驗(yàn)中將ROI內(nèi)信號(hào)強(qiáng)度平均值M(ROI)的3倍作為篩選X射線噪聲點(diǎn)的閾值，以便能有效的篩選出X射線噪聲點(diǎn)。

(2)對(duì)噪聲點(diǎn)進(jìn)行賦值(中值濾波算法)。由于X射線噪聲點(diǎn)多為孤立噪聲，篩選出噪聲點(diǎn)后，可以3×3為模板，使用中值濾波算法為X射線噪聲點(diǎn)賦值。

1.3 雙波長(zhǎng)減影算法

非目標(biāo)熒光信號(hào)屬于確定性信號(hào)而非隨機(jī)噪聲，在不同熒光濾光片作用下產(chǎn)生的熒光信號(hào)其性質(zhì)相同，只是在強(qiáng)度上存在一個(gè)比例關(guān)系。因此，可以通過(guò)雙波長(zhǎng)減影算法來(lái)消除非目標(biāo)熒光團(tuán)信號(hào)，其雙波長(zhǎng)圖像減影算法為公式3：

式中SX為待求解的目標(biāo)熒光信號(hào)，k為減影系數(shù)，Sλ1和Sλ2為分別在λ1和λ2波長(zhǎng)采集的圖像，包含目標(biāo)熒光信號(hào)和非目標(biāo)熒光信號(hào)。

減影系數(shù)k是一個(gè)與實(shí)驗(yàn)設(shè)置(如濾光片的中心波長(zhǎng)、帶寬、透過(guò)率以及CCD相機(jī)的光譜靈敏度等)有關(guān)的常數(shù)，減影系數(shù)k可以根據(jù)不同波段采集的熒光圖中相同區(qū)域X內(nèi)的信號(hào)的相對(duì)強(qiáng)度求解，其最小二乘問(wèn)題的求解為公式4：

2 仿真實(shí)驗(yàn)

使用仿真實(shí)驗(yàn)來(lái)對(duì)算法的可行性進(jìn)行驗(yàn)證。

2.1 仿真圖像構(gòu)建

構(gòu)建像素500×500大小的仿真圖像，區(qū)域1像素150×150模擬目標(biāo)熒光團(tuán)，區(qū)域2像素150×150模擬非目標(biāo)熒光團(tuán)，在不同熒光波長(zhǎng)λ1，λ2下，區(qū)域1分別賦值14000和12000，區(qū)域2分別賦值3000和24000。區(qū)域1及2為ROI，ROI以外區(qū)域?yàn)楸尘皡^(qū)域。整幅圖的環(huán)境噪聲強(qiáng)度賦值200，隨機(jī)選取25個(gè)像素點(diǎn)賦值為63000，模擬高強(qiáng)度X射線噪聲。為了更接近真實(shí)，整幅圖像加入20 db高斯白噪聲，如圖1所示。

圖1 構(gòu)建的λ1和λ2仿真圖像

2.2 仿真圖像預(yù)處理

求取背景區(qū)域的像素點(diǎn)均值，作為環(huán)境噪聲閾值m。在本次仿真實(shí)驗(yàn)中，得到m=205，然后采用硬閾值算法進(jìn)行環(huán)境噪聲去除。

改進(jìn)型中值濾波算法去除X射線噪聲的流程為：①使用閾值法對(duì)X射線強(qiáng)噪聲點(diǎn)進(jìn)行篩選，根據(jù)2.2部分的理論推導(dǎo)，采用ROI內(nèi)數(shù)值平均值的3倍作為閾值P，大于P的元素都認(rèn)為是強(qiáng)噪聲點(diǎn)。本例中可以求得P=31733，共篩選出25個(gè)強(qiáng)噪聲點(diǎn)。應(yīng)用中值濾波算法，對(duì)篩選出的噪聲點(diǎn)進(jìn)行3×3模板的算法處理。

經(jīng)過(guò)圖像預(yù)處理之后，背景區(qū)域的環(huán)境噪聲和隨機(jī)分布的X射線高強(qiáng)度噪聲點(diǎn)被有效濾除，且很好地保留圖像的邊緣特性。預(yù)處理前后的2D及3D對(duì)比如圖2所示。

圖2 λ1仿真圖像預(yù)處理前后對(duì)比圖

2.3 仿真圖像減影處理

熒光波長(zhǎng)λ1、λ2圖像分別經(jīng)過(guò)預(yù)處理之后，本研究進(jìn)一步開(kāi)始減影處理。在兩幅圖的非目標(biāo)熒光區(qū)域2內(nèi)選取5個(gè)像素10×10的模板求取減影系數(shù)k。5個(gè)模板的k值分別為0.508、0.504、0.498、0.496及0.504，取平均值0.5020為本次實(shí)驗(yàn)的減影系數(shù)，如圖3所示。

圖3 計(jì)算減影系數(shù)的模板分布圖

由公式(3)用λ1圖像減去k倍的λ2圖像，得到減影后結(jié)果。λ1圖像減影前后的2D及3D對(duì)比如圖4所示。從圖中可以直觀地看出，去除了非目標(biāo)熒光信號(hào)的干擾，同時(shí)有效地保留了目標(biāo)熒光信號(hào)。

2.4 效果評(píng)估

本研究使用目標(biāo)與背景信號(hào)強(qiáng)度之比(target-tobackground ratio，TBR)來(lái)定量評(píng)估所提圖像處理方法的性能[12-13]。目標(biāo)信號(hào)是指目標(biāo)區(qū)域內(nèi)的信號(hào)，背景信號(hào)是指目標(biāo)區(qū)域外的信號(hào)，因此TBR即目標(biāo)區(qū)域內(nèi)的信號(hào)強(qiáng)度與目標(biāo)區(qū)域外的信號(hào)強(qiáng)度之比，其TBR的計(jì)算為公式5：

圖4 λ1仿真圖像減影前后對(duì)比圖

式中T表示目標(biāo)區(qū)域內(nèi)各像素信號(hào)強(qiáng)度總和，B表示目標(biāo)區(qū)域外各像素信號(hào)強(qiáng)度總和。

本次仿真實(shí)驗(yàn)中XLI圖像處理前后的TBR對(duì)照見(jiàn)表1。

表1 XLI圖像處理前后TBR對(duì)照

結(jié)果表明，應(yīng)用所提方法進(jìn)行圖像處理之后，圖像TBR較處理前提高約38倍，表明了本次仿真實(shí)驗(yàn)中，經(jīng)過(guò)硬閾值算法、中值濾波算法和雙波長(zhǎng)減影算法處理后，圖像質(zhì)量得到了很好的提升。

3 仿體實(shí)驗(yàn)

為了進(jìn)一步評(píng)估算法的性能，本研究開(kāi)展仿體實(shí)驗(yàn)研究，其處理方法與仿真實(shí)驗(yàn)部分類(lèi)似。仿體實(shí)驗(yàn)所使用的圖像采集系統(tǒng)為XLI系統(tǒng)[14-15]。

3.1 仿體設(shè)置

成像物體采用的是高40 mm，直徑為30 mm的圓柱體，如圖5(a)所示。在1號(hào)試管內(nèi)加入Y2O3：Eu作為目標(biāo)熒光團(tuán)，在2號(hào)試管內(nèi)加入Gd2O2S：Tb作為非目標(biāo)熒光團(tuán)。試管直徑約4 mm，2個(gè)試管邊邊距約為7 mm。感興趣區(qū)域?yàn)槌上裎矬w中空部分(用紅色虛線框標(biāo)示)，整幅圖的ROI以外部分為背景區(qū)域。

使用50 kv，1 mA的X射線激發(fā)探針，目標(biāo)熒光團(tuán)發(fā)紅光(中心波長(zhǎng)約為600 nm)，非目標(biāo)熒光團(tuán)發(fā)綠光(中心波長(zhǎng)約為540 nm)，以(600±37)nm和(520±44) nm的帶通濾光片進(jìn)行濾光，通過(guò)XLI系統(tǒng)可以分別獲得2個(gè)波段的XLI圖像。圖5(b)為以600 nm為中心波長(zhǎng)的XLI圖像，從圖中可以看出存在高強(qiáng)度的X射線噪聲和非目標(biāo)熒光信號(hào)(Gd2O2S：Tb)的干擾，降低了目標(biāo)熒光信號(hào)(Y2O3：Eu)的可信度和可辨識(shí)度。

圖5 仿體設(shè)置及中心波長(zhǎng)600 nm的XLI圖像

3.2 仿體圖像預(yù)處理

為了更直觀的觀察圖像及處理結(jié)果，對(duì)圖像采用偽彩顯示。以波長(zhǎng)600 nm下XLI圖像為例來(lái)闡述仿體圖像的預(yù)處理過(guò)程。

(1)以背景區(qū)域環(huán)境噪聲強(qiáng)度平均值(m=363)為閾值進(jìn)行硬閾值算法處理，去除環(huán)境噪聲。

(2)去除X射線噪聲，求得感興趣區(qū)域內(nèi)信號(hào)強(qiáng)度平均值的3倍(P=25843)為篩選X射線噪聲點(diǎn)的閾值，共篩選出36個(gè)高強(qiáng)度噪聲點(diǎn)。對(duì)篩選出的噪聲點(diǎn)進(jìn)行3×3模板的中值濾波處理。預(yù)處理前后的2D及3D對(duì)比如圖6所示，可以直觀地看出，預(yù)處理去除了環(huán)境噪聲及高強(qiáng)度X射線噪聲點(diǎn)。

圖6 600 nm波段的XLI圖像的預(yù)處理前后對(duì)比圖

3.3 仿體圖像減影處理

本次實(shí)驗(yàn)的非目標(biāo)熒光主要是指2號(hào)試管的熒光信號(hào)。考慮到非目標(biāo)熒光團(tuán)區(qū)域內(nèi)不同位置熒光信號(hào)的相對(duì)強(qiáng)度可能會(huì)略有不同，本研究在非目標(biāo)區(qū)域內(nèi)隨機(jī)選擇5個(gè)像素10×10區(qū)域作為模板進(jìn)行減影系數(shù)求解，模板具體位置用黑色實(shí)線方框標(biāo)示，如圖7所示。

圖7 計(jì)算減影系數(shù)的模板分布圖

由公式(4)求得的5個(gè)模板的k值分別為0.56、0.59、0.57、0.59以及0.59，取其平均值0.58為本次實(shí)驗(yàn)的減影系數(shù)進(jìn)行減影處理。

以600 nm下XLI圖像為例，減影前后的2D及3D對(duì)比如圖8所示。從圖中可以看出非目標(biāo)熒光被濾除，而目標(biāo)熒光得以保留。

圖8 600 nm采集圖像的減影前后對(duì)比圖

3.4 效果評(píng)估

波長(zhǎng)600 nm下的XLI圖像處理前后的對(duì)比如圖9所示。從圖中可以直觀地看出，經(jīng)過(guò)預(yù)處理以及雙波長(zhǎng)減影算法處理后，環(huán)境噪聲、X射線噪聲以及非目標(biāo)熒光信號(hào)均被有效濾除。

圖9 600 nm下XLI圖像處理前后對(duì)比圖

進(jìn)一步計(jì)算TBR來(lái)定量考察所提方法的性能，圖像處理之后的TBR比處理前提高約4倍，見(jiàn)表2。

表2 圖像處理前后TBR對(duì)比

4 結(jié)論

XLI圖像容易受到環(huán)境噪聲、泄露的X射線以及非目標(biāo)熒光等背景噪聲干擾，造成圖像信噪比較低，從而導(dǎo)致成像目標(biāo)的可辨識(shí)度和可信度不高。在本研究中實(shí)現(xiàn)了一種基于雙波長(zhǎng)圖像減影的目標(biāo)熒光信號(hào)提取方法，首先根據(jù)不同噪聲的特點(diǎn)采用硬閾值算法及改進(jìn)型中值濾波算法對(duì)圖像進(jìn)行預(yù)處理，去除環(huán)境噪聲和X射線強(qiáng)噪聲干擾，然后基于目標(biāo)熒光探針和非目標(biāo)熒光探針發(fā)光光譜的不同，結(jié)合最小二乘法確定減影系數(shù)，并進(jìn)行雙波長(zhǎng)圖像減影處理，去除非目標(biāo)熒光信號(hào)，提取目標(biāo)熒光信號(hào)。本研究通過(guò)仿真實(shí)驗(yàn)和仿體實(shí)驗(yàn)對(duì)所提方法進(jìn)行了有效性驗(yàn)證，效果評(píng)估結(jié)果表明，所提方法可以有效濾除圖像中的背景噪聲，提高圖像的TBR，從而提高了成像質(zhì)量，有利于成像目標(biāo)的識(shí)別。

[1]Pratx G，Carpenter CM，Sun C，et al.X-ray luminescence computed tomography via selective excitation：A feasibility study[J].IEEE Transactions on Medical Imaging，2010，29(12)：1992-1999.

[2]Pratx G，Carpenter CM，Sun C，et al.Tomographic molecular imaging of X-ray-excitable nanoparticles[J].Optics Letter，2010，35(20)：3345-3347.

[3]Pratx G，Carpenter CM，Sun C，et al.Hybrid X-ray/optical luminescence imaging：characterization of experimental conditions[J].Medical Physics，2010，37(8)：4011-4018.

[4]Weissleder R，Pittet MJ.Imaging in the era of molecular oncology[J].Nature，2008，452(7187)：580-589.

[5]Willmann JK，Bruggen NV，Dinkelborg LM，et al. Molecular imaging in drug development[J]. Nature Reviews Drug Discovery，2008，7(7)：591-607.

[6]Lv Y，Tian J，Cong W，et al.A multilevel adaptive finite element algorithm for bioluminescence tomography[J].Optics Express，2006，14(18)：8211-8223.

[7]Naczynski DJ，Sun C，Türkcan S，et al.X-ray-induced shortwave infrared biomedical imaging using rare-earth nanoprobes[J].Nano Letters，2015，15(1)：96-102.

[8]Liu F，Liu X，Zhang B，et al.Extraction of target fluorescence signal from in vivo background signal using image subtraction algorithm[J]. International Journal of Automation and computing，2012，9(3)：232-236.

[9]Deliolanis NC，Wurdinger T，Pike L，et al.In vivo tomographic imaging of red-shifted fluorescent proteins[J].Biomedical Optics Express，2011，2(4)：887-900.

[10]Chen H，Moore T，Qi B，et al.Monitoring pH-triggered drug release from radioluminescent nanocapsules with X-ray excited optical luminescence[J].ACS nano，2013，7(2)：1178-1187.

[11]Silva AAD，Cebim MA and Davolos MR. Excitation mechanisms and effects of dopant concentration in Gd2O2S：Tb3+phosphor[J]. Journal of Luminescence，2008，128(7)：1165-1168.

[12]Wang P，Meng Z，Tan J，et al.An improved method for measurement of target-to-background ratio in assessing mediastinal lesions on18F-FDG coincidence SPECT/CT imaging[J]. Nuclear Medicine Communications，2010，31(5)：398-404.

[13]Yang C，Hou VW，Girard EJ，et al.Targetto-background enhancement in multispectral endoscopy with background autofluorescence mitigation for quantitative molecular imaging[J]. Journal of Biomedical Optics，2014，19(7)：76014.

[14]Liu X，Liao QM，Wang HK.Fast X-Ray Luminescence Computed Tomography Imaging[J]. IEEE Transactions on Biomedical Engineering，2014，61(6)：1621-1627.

[15]Liu X，Liao QM，Wang HK，et al.Excitationresolved cone-beam x-ray luminescence tomography[J].Journal of Biomedical Optics，2015，20(7)：070501.

Research of X-ray luminescence imaging based on dual wavelength subtraction algorithm/

NING Qian, XU Huan, LIU Yang, et al//
China Medical Equipment,2017,14(9):6-10.

Objective: To discuss the removal effectiveness of dual wavelength subtraction algorithm for background noises, including environmental noises, X-ray noises and nontarget luminescence signals and so on, which existed in X-ray luminescence imaging(XLI). Methods: The XLI processing method based on dual-wavelength subtraction algorithm was developed, and it was verified through numerical simulation and phantom experiment. The original image was preprocessed through using a hard threshold algorithm and improved median filter algorithm, and then the dual-wavelength subtraction algorithm was further used to conduct image. Finally, the performance of the proposed results were verified and assessed. Results: The target-to-background ratio(TBR) of XLI was enhanced to approximately 38 times in numerical simulation and it was enhanced to approximately 4 times in phantom experiment. The experiments showed that this method could increased TBR of XLT imaging and the imaging quality were improved. Conclusion: The proposed image processing method which based on dual wavelength subtraction algorithm can efficiently remove the interference of background noises and effectively extract the luminescence signals of target.

X-ray luminescence imaging; Optical molecular imaging; Background noise; Dual-wavelength subtraction; Target-to-background ratio

School of Biomedical Engineering, Fourth Military Medical University, Xi’an 710032, China.

1672-8270(2017)09-0006-05

R812

A

10.3969/J.ISSN.1672-8270.2017.09.002

寧騫,男,(1995- ),本科學(xué)歷。第四軍醫(yī)大學(xué)生物醫(yī)學(xué)工程學(xué)院，研究方向：醫(yī)學(xué)成像與圖像處理。

2017-06-22

陜西省社會(huì)發(fā)展科技攻關(guān)項(xiàng)目(2016SF-044)“多能級(jí)X射線激發(fā)的多光譜熒光斷層成像技術(shù)研究”

①第四軍醫(yī)大學(xué)生物醫(yī)學(xué)工程學(xué)院 陜西 西安 710032

②中央軍委后勤保障部衛(wèi)生局藥品儀器檢驗(yàn)所 北京 100071

*通訊作者：que_sera@yeah.net