白敏茹,黃孝龍,顧廣澤,趙雪瑩
(湖南大學(xué) 數(shù)學(xué)與計量經(jīng)濟學(xué)院,湖南 長沙 410082)
?
基于張量秩校正的圖像恢復(fù)方法
白敏茹*,黃孝龍,顧廣澤,趙雪瑩
(湖南大學(xué) 數(shù)學(xué)與計量經(jīng)濟學(xué)院,湖南 長沙 410082)
針對醫(yī)學(xué)圖像和視頻圖像的恢復(fù)問題,基于張量表示,研究有限樣本下的低秩張量數(shù)據(jù)恢復(fù)問題,在張量奇異值分解(t-SVD)理論的基礎(chǔ)上,提出了張量秩校正模型和兩階段張量秩校正方法,第一階段是用張量核范數(shù)最小化模型求得預(yù)估解,第二階段,根據(jù)預(yù)估解,求解張量秩校正模型,獲得更高精度的解.構(gòu)建了求解張量秩校正模型和張量核范數(shù)最小化模型的張量近似點算法,使得可以在實數(shù)域上對張量直接進行計算,并且從理論上證明了該算法的收斂性.通過對醫(yī)學(xué)圖像和視頻圖像的數(shù)值仿真實驗,驗證了本文所提出模型和方法的有效性,實驗結(jié)果顯示,張量秩校正模型和方法能夠取得更高的恢復(fù)精度.
圖像恢復(fù);張量奇異值分解;張量秩校正;張量近似點算法
隨著電子技術(shù)和成像技術(shù)的發(fā)展,從醫(yī)學(xué)圖像到遙感圖像,從導(dǎo)彈精確制導(dǎo),到人臉識別及指紋識別再到具有視覺功能的智能機器人,人類活動的方方面面都會產(chǎn)生或涉及到大量的高維圖像.高維圖像已經(jīng)成為一種重要的多媒體形式,廣泛存在于人們的日常生活中.圖像在形成,傳輸和記錄的過程中受多種因素的影響,圖像的質(zhì)量會有所下降,典型表現(xiàn)為色彩模糊和有噪聲干擾等.這一降質(zhì)的過程被稱為圖像的退化.圖像恢復(fù)的目的就是盡可能地恢復(fù)退化了的高維圖像的本來面目.
傳統(tǒng)的圖像處理方法是基于向量和矩陣的表示形式,往往破壞了這些數(shù)據(jù)的原始空間結(jié)構(gòu),在分析過程中不能夠很好地刻畫這些數(shù)據(jù)的本質(zhì)和充分挖掘其內(nèi)部特性.張量作為向量和矩陣表示的高階推廣,能夠更好地表達高階數(shù)據(jù)復(fù)雜的本質(zhì)結(jié)構(gòu),已被廣泛應(yīng)用于計算機視覺與圖像、人臉識別、醫(yī)學(xué)圖像和統(tǒng)計信號處理等研究領(lǐng)域中[1-6].
高維圖像數(shù)據(jù)往往具有低維屬性,張量完備化問題就是利用張量數(shù)據(jù)的低秩結(jié)構(gòu),是一種在有限樣本或測量數(shù)據(jù)下最小化張量的秩的優(yōu)化問題.最小化張量的秩是NP難問題,通常的處理方法有:1)將張量轉(zhuǎn)化成矩陣,然后求解矩陣完備化問題[7];2)用特殊的張量分解方法來分解張量,如CANDECOMP/PARA-FAC(CP)分解,Tucker分解等方法.
由于矩陣的核范數(shù)是矩陣秩的緊的凸逼近,因此對矩陣完備化問題的求解一般是將其轉(zhuǎn)化為矩陣核范數(shù)最小化問題求解.對矩陣核范數(shù)最小化問題的求解有近似點算法(PPA)[8],交替方向方法(ADM),加速近似梯度方法(APG)[9].雖然低秩矩陣完備化問題得到很好發(fā)展,但張量完備化問題研究還很不完善.不同于矩陣秩只有一種定義,張量秩有多種定義.傳統(tǒng)上主要有兩種張量秩的定義,CP秩和Tucker秩,它們分別是基于CP分解和 Tucker分解的.將張量展開成矩陣,利用展開矩陣性質(zhì)近似逼近張量的秩,是常用的處理方法.例如:Gandy[2]等用各片分別展開矩陣的核范數(shù)的和作為張量秩的近似逼近;Liu[5]等進一步將各片分別展開矩陣的核范數(shù)通過加權(quán)來近似張量的秩,并提出了HaLTRC算法求解該松弛模型(TSN).然而這兩種逼近方法并不是張量秩函數(shù)的最緊的凸逼近[7].
Kilmer等[10]基于快速傅里葉變換可以將塊循環(huán)矩陣對角化的思想,提出了張量奇異值分解(T-SVD)方法,使得張量可以在傅里葉變換下實現(xiàn)快速分解.基于T-SVD, Semerci等[6]提出張量核范數(shù)概念,對于3階張量,利用張量核范數(shù)近似逼近張量的秩,建立了張量核范數(shù)最小化模型(TNN),構(gòu)建了交替方向方法(ADMM)求解該模型,并應(yīng)用于多線性數(shù)據(jù)的圖像壓縮和恢復(fù),通過對比,TNN逼近比TSN逼近效果更好.但是該文沒有給出ADMM方法的收斂性結(jié)果,文中的ADMM算法一部分在實數(shù)域上計算,一部分在復(fù)數(shù)域上計算.與以往模型不一樣,TNN模型的目標(biāo)變量是定義在復(fù)數(shù)域即傅里葉域內(nèi)的矩陣,約束變量是定義在實數(shù)域的.因此,根據(jù)這個問題的特點,設(shè)計更加有效的具有收斂性的優(yōu)化算法,是亟需解決的一個問題.另外,文獻[11]指出,矩陣核范數(shù)在某些情況下不是矩陣秩的最緊凸逼近,如對角元素被高度樣本化,則矩陣核范數(shù)最小化模型求解低秩恢復(fù)問題的能力就會高度弱化,而矩陣核范數(shù)是張量核范數(shù)(TNN)的二階形式.本文針對以上兩個問題開展研究,主要貢獻有兩個:一是提出了張量秩校正模型(CRTNN)和兩階段張量秩校正方法,二是構(gòu)建了張量近似點算法,用于求解CRTNN模型和TNN模型,從理論上證明了該算法的收斂性.仿真實驗驗證了本文所提出模型和方法的有效性.結(jié)果顯示,在醫(yī)學(xué)圖像以及視頻圖像的恢復(fù)問題中,張量秩校正方法能夠取得更高的恢復(fù)精度.
張量即為多維數(shù)組,其元素所在位置需要3個或3個以上的變量來表示,可以記為A=an1n2…nN,其中A∈Rn1×n2×…×nN,這里ni,i=1,2, …,nN,稱為維數(shù),N為階數(shù).特別地,向量為一維張量,矩陣為二維張量.張量A的片是只有兩個指標(biāo)沒有固定的矩陣形式,例如3階張量A∈Rn1×n2×n3的前片表示為A(:,:,i),i=1,2,3.針對彩色圖像的數(shù)據(jù)結(jié)構(gòu),本文主要討論3階張量.本節(jié)簡要介紹與張量奇異值分解有關(guān)的基本知識,更多詳細的介紹見文獻[10].
首先介紹由張量的前片構(gòu)成的塊循環(huán)矩陣,對張量A∈Rn1×n2×n3,其前片為n1×n2的矩陣A1,A2,A3其中Ai=A(:,:,i),則有:
(1)
而快速傅里葉變換(FFT)可以將塊循環(huán)矩陣轉(zhuǎn)化成塊對角矩陣,有如下形式:
(2)
(3)
(4)
取MatVec變換作用于張量A∈Rn1×n2×n3的每一個前片,則MatVec(A)將張量A∈Rn1×n2×n3作用成n1n3×n2的矩陣:
(5)
將MatVec(A)返回張量則用fold變換:
fold(MatVec(A))=A.
(6)
定義1.1[10]張量A∈Rn1×n2×n3和B∈Rn2×l×n3,則A與B的張量積定義如下:
C:=A*B=fold(circ(A)·MatVec(B)).
(7)
其中C∈Rn2×l×n3.
接下來介紹單位張量,張量的轉(zhuǎn)置,張量奇異值分解(t-SVD).
定義1.2[10]張量In×n×l∈Rn×n×l為單位張量當(dāng)它的第一個前片為n×n的單位矩陣,剩余兩前片的元素均為0.
定義1.3[10]張量A∈Rn1×n2×n3,則張量A的張量轉(zhuǎn)置AT∈Rn2×n1×n3,是通過保持第1個前片位置不變,第2個前片和第3個前片位置互換,并且分別對張量A每一個前片做轉(zhuǎn)置得來.
定義1.4[10]張量A∈Rn×n×l為f-對角化張量當(dāng)且僅當(dāng)其每一個前片矩陣均為對角矩陣.
定義1.5[10]張量A∈Rn×n×l為正交張量當(dāng)且僅當(dāng)它滿足:
AT*A=A*AT=In×n×l.
(8)
定理 1.6[10](張量奇異值分解(t-SVD))張量A∈Rn1×n2×n3,則A可以分解為:
A=U*S*VT.
(9)
其中U和V為正交張量,其中U∈Rn1×n1×n3,V∈Rn2×n2×n3,S∈Rn1×n2×n3為f-對角化張量.
最后介紹三階張量的多線性秩和張量核范數(shù)(TNN)以及它們之間的聯(lián)系.
在一定條件下,矩陣核范數(shù)是矩陣的秩的凸松弛,同理,張量核范數(shù)(TNN)是張量多線性秩的凸松弛[13].
在圖像采集過程中,由于各種原因,可能會出現(xiàn)圖像數(shù)據(jù)損失的情況,即圖像序列組成的張量X中有部分元素的值缺失.張量數(shù)據(jù)恢復(fù)問題即為張量完備化問題,就是利用張量數(shù)據(jù)的低秩結(jié)構(gòu),在有限樣本或測量數(shù)據(jù)下最小化張量的秩的優(yōu)化問題.Liu[5]等將張量按照不同方向分別展開成矩陣的核范數(shù)通過加權(quán)來近似張量的秩,建立了張量完備化的如下凸松弛模型(TSN):
(10)
其中X(i)是張量X 的i模矩陣,并提出了HaLTRC算法求解該TSN模型.
Semerci等[6]則基于張量核范數(shù)TNN,提出了張量核范數(shù)最小化模型(TNN):
(11)
這里張量X,M∈Rn1×n2×n3,M在集合Ω里的元素是給定的,不在Ω里的元素則是0,即:
(12)
注意到矩陣核范數(shù)在某些情況下不是矩陣秩的最佳凸逼近,如對角元素被高度樣本化,則矩陣核范數(shù)最小化模型求解低秩恢復(fù)問題的能力會高度弱化[11].為了獲得更高精度解,針對張量核范數(shù)最小化模型(TNN),本文提出了一個張量秩校正模型(CRTNN):
s.t.XΩ=MΩ.
(13)
案例3:女,59歲,于2017年5月前因乳房疼痛不適,在外院檢查發(fā)現(xiàn)乳腺癌,病理活檢診斷為非特殊性浸潤癌,ER(-),PR(-),HER2(-),行全身檢查,準(zhǔn)備手術(shù)治療,檢查發(fā)現(xiàn)顱內(nèi)多發(fā)占位,考慮乳癌腦轉(zhuǎn)移,預(yù)后不佳,不建議乳癌手術(shù)治療,給予相關(guān)保守治療。后經(jīng)幾次放療,效果不佳,患者漸漸意識狀況下降,經(jīng)綜合考慮,患者家人要求免疫治療。2018年3月在我院進行免疫治療,經(jīng)過3個多月治療患者腦內(nèi)腫瘤未見明顯增大,患者各種神經(jīng)系統(tǒng)癥狀均予以改善,KPS評分值提高。
F(D)=UDiag(f(σ(D)))VT.
σ(D)為矩陣D的奇異值,對稱函數(shù)f:Rn→R定義如下:
對于τ>0,ε>0,標(biāo)量函數(shù)φ:R→R表示如下:
其中sgn(t)為符號函數(shù),定義如下:
(14)
針對張量秩校正模型,下面給出兩階段張量秩校正方法:
3.1 算法求解
在本節(jié)中,應(yīng)用近似點算法(PPA)求解式(11)和式(13).
問題(13)通過轉(zhuǎn)化為式(14),應(yīng)用近似點算法(PPA)求解.
式(14)的拉格朗日函數(shù)為:
則近似點算法(PPA)求解:
(15)
其中γ∈(0,2),r·s>1.
(16)
(17)
(18)
(19)
由于式(15)是在傅里葉域里求解式 (14)矩陣核范數(shù)最小化問題.然而原問題式(13)是在實數(shù)域里的結(jié)果,并且變量都是以張量形式計算.考慮到逆傅里葉變換和張量奇異值分解(t-SVD)的思想,在每一步迭代中應(yīng)用如下步驟:
通過上述變換,則可以把矩陣迭代過程轉(zhuǎn)換為張量迭代過程:
(20)
將式(20)里的矩陣形式轉(zhuǎn)化為張量形式,并還回實數(shù)域:
(21)
(22)
綜上,得出求解問題(12)的近似點算法(PPA)如下:
算法3.1
任取γ∈(0,2),及r·s>1,對給定(Xk,Yk) ,
.(23)
2)松弛步:新的迭代步,
(24)
3.2 收斂性分析
定理3.1 設(shè){Xk,Yk}為由算法3.1產(chǎn)生的迭代序列,若γ∈(0,2),及rs>1,則序列{Xk,Yk}收斂到問題(13)的鞍點.
證明 對于矩陣求解問題,由He[8]定理3.7分析,當(dāng)k→時,易知:
(25)
(26)
由于快速傅里葉變換和逆快速傅里葉變換均為連續(xù)有界算子,故有:
(27)
取MatVec(Xk+1),MatVec(X*),MatVec(Yk+1)以及MatVec(Y*)分別表示circ(Xk+1),circ(X*),circ(Yk+1),和circ(Y*)的第一列塊循環(huán)矩陣,并作fold變換.則得到:
(28)
X*和Y*為(13)的鞍點.
證畢.
本文針對醫(yī)學(xué)圖像和視頻圖像的恢復(fù)問題,分別對張量核范數(shù)加權(quán)和模型(TSN模型)式(10)、張量核范數(shù)最小化模型(TNN模型)式(11)以及張量秩校正模型(CRTNN模型)式(13)進行仿真實驗,TSN模型采用HaLRTC方法求解,TNN模型和CRTNN模型采用近似點算法(PPA)求解,并給出仿真結(jié)果,所有結(jié)果都是在Core i5 的CPU及4G內(nèi)存的Windows 7系統(tǒng)下的ASUS筆記本中運行MATLAB R2012b計算得出.圖像恢復(fù)的數(shù)值評價指標(biāo)通常由相對誤差和峰值信噪比(PSNR)計算,相對誤差計算公式:
(29)
式中:M為實值張量;X為預(yù)估張量,峰值信噪比計算公式:
(30)
式中:n1,n2,n3為張量M∈Rn1×n2×n3的維數(shù),同時終止條件為:
(31)
tol為終止參數(shù),取tol=10-3,主要是小于這個值之后,變化特別微小.
文中選取的圖像為大小415×477×3的醫(yī)學(xué)圖像,視頻圖像為大小112×160×3的視頻的其中一幀,進行仿真實驗,并比較TSN模型,TNN模型,CRTNN模型的恢復(fù)效果.
圖1為醫(yī)學(xué)圖像和視頻圖像原始圖像.圖2,圖3分別為醫(yī)學(xué)圖像和視頻圖像在樣本率為20%(即有效信息只有20%)的情況時用TSN模型,TNN模型,CRTNN模型視覺恢復(fù)效果對比,從圖2,圖3的PSNR值對比和視覺恢復(fù)效果對比中,可以發(fā)現(xiàn)本文提出的CRTNN模型能得到更好的恢復(fù)效果.
圖4分別為醫(yī)學(xué)圖像和視頻圖像在TSN模型,TNN模型,CRTNN模型下對不同樣本率得到的相對誤差曲線對比.從中可以明顯看出:本文提出的張量秩校正方法對不同的樣本率得到的恢復(fù)圖像的相對誤差曲線都是最低的,表明本文提出的CRTNN模型能夠取得更高精度的恢復(fù)效果.
圖1 原始圖像
圖2 樣本率為20%的醫(yī)學(xué)圖像及其分析在TSN模型,TNN模型和CRTNN模型下的恢復(fù)圖
Fig.2 Recovery results on a medical image with 20% sample ratio by models TSN,TNN and CRTNN, respectively
圖3 樣本率為20%的視頻圖像及其分別在TSN模型,TNN模型和CRTNN模型下的恢復(fù)圖
Fig.3 Recovery results on a video image with 20% sample ratio by models TSN, TNN and CRTNN, respectively
(a)醫(yī)學(xué)圖像結(jié)果
(b)視頻圖像結(jié)果
針對高維圖像恢復(fù)問題,本文提出了張量秩校正模型和兩階段張量秩校正方法,并提出了求解張量秩校正模型的張量近似點算法,從理論上分析了該算法的收斂性.仿真結(jié)果驗證了本文所提出模型和方法的有效性,結(jié)果表明,張量秩校正方法模型能夠取得更高的恢復(fù)精度.能否將該模型和算法推廣到四階及以上的圖像恢復(fù)問題?這個問題值得進一步研究.
[1] ELY G, AERON S, MILLER E L. Exploiting structural complexity for robust and rapid hyper spectral imaging [C]//Proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2013:2193-2197.
[2] GANDY S, RECHT B, YAMADA I. Tensor completion and low-n-rank tensor recovery via convex optimization [J]. Inverse Problems, 2011, 27(2): 025010.
[3] HAO N H, KILMER M E, BRAMAN K, et al.Facialrecognitionwithtensor-tensordecompositions[J].SIAMJournalonImagingSciences, 2013, 6(1): 437-463.
[4]KILMERM,BRAMANK,HAON,et al.Third-ordertensorsasoperatorsonmatrices:Atheoreticalandcomputationalframeworkwithapplicationsinimaging[J].SIAMJournalonMatrixAnalysisandApplications, 2013, 34(1):148-172.
[5]LIUJ,MUSIALSKIP,WONKAP, et al.Tensorcompletionforestimatingmissingvaluesinvisualdata[J].IEEETransactionsonPatternAnalysisandMachineIntelligence, 2013, 35(1): 208-220.
[6]SEMERCIO,HAON,KILMERME,et al.Tensor-basedformulationandnuclearnormregularizationformultienergycomputedtomography[J].IEEETransactiononImageProcessing, 2014, 23(4): 1678-1693.
[7]MUC,HUANGB,WRIGHTJ,et al.Squaredeal:Lowerboundsandimprovedrelaxationsfortensorrecovery[C]//Proceedingsofthe31stInternationalConferenceonMachineLearning(ICML-14), 2014, 32(1): 73-81.
[8]HEBS,YUANXM,ZHANGWX.Acustomizedproximalpointalgorithmforconvexminimizationwithlinearconstraints[J].ComputationalOptimizationandApplications, 2013, 56(3): 559-572.
[9]TOHKC,YUNS.Anacceleratedproximalgradientalgorithmfornuclearnormregularizedleastsquaresproblems[J].PacificJournalofOptimization, 2010, 6(3): 615-640.
[10]KILMERME,MARTINCD.Factorizationstrategiesforthird-ordertensors[J].LinearAlgebraanditsApplications, 2011, 435(3):641-658.
[11]MIAOW,PANS,SUND.Arank-correctedprocedureformatrixcompletionwithfixedbasiscoefficients[J].Math.Programming,2016,159(1):289-338.
[12]ZHANGZ,ELYG,AERONS,et al.Novelmethodsformultilineardatacompletionandde-noisingbasedontensor-SVD[C]//InProceedingsoftheIEEEConferenceonComputerVisionandPatternRecognition(CVPR), 2014, 3842-3849.
[13]CAIJF,CANDESEJ,SHENZ.Asingularvaluethresholdingalgorithmformatrixcompletion[J].SIAMJournalonOptimization, 2010, 20(4): 1956-1982.
Tensor Rank Corrected Procedure for Image Restoration
BAI Min-ru?,HUANG Xiao-long,GU Guang-ze,ZHAO Xue-ying
(College of Mathematics and Econometrics, Hunan Univ, Changsha, Hunan 410082, China )
Tensor-based restoration of medical images and video images was studied with limited samples. On the basis of the theory of tensor singular value decomposition (t-SVD), a tensor rank-correction model (CRTNN) was proposed to correct the tensor nuclear norm minimization model (TNN). A two-stage rank correction method is given as follows: the first stage is used to generate a pre-estimator by solving the TNN model, and the second stage is to solve the CRTNN model to generate a high-accuracy recovery by the pre-estimator. A tensor proximal point algorithm was proposed to solve the CRTNN model and the TNN model, making it possible to calculate tensor directly in the real field. The convergence of the algorithm was proved in theory. Numerical experiments of medical images and video images verify the efficiency of the proposed model and method. The experiment results show that tensor rank-correction model and method can achieve higher-accuracy recovery.
image restoration;t-SVD; tensor rank-correction model; tensor proximal point algorithm
1674-2974(2016)10-0148-07
2016-01-17
國家自然科學(xué)基金資助項目(11571098), National Natural Science Foundation of China(11571098) ; 湖南省高校創(chuàng)新平臺開放基金資助項目(14K018)
白敏茹(1968-),女,江西宜春人,湖南大學(xué)博士生導(dǎo)師, 副教授
?通訊聯(lián)系人,E-mail:minru-bai@163.com
TP751
A