• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    總變差約束的數(shù)據(jù)分離最小圖像重建模型及其Chambolle-Pock求解算法?

    2018-11-03 04:32:18喬志偉
    物理學報 2018年19期
    關鍵詞:模體范數(shù)高精度

    喬志偉

    (山西大學計算機與信息技術學院,太原 030006)

    (2018年4月27日收到;2018年7月2日收到修改稿)

    1 引 言

    在計算機斷層成像(computed tomography,CT)圖像重建中,以濾波反投影(filtered backprojection,FBP)法為代表的解析重建方法在商用CT中占據(jù)主流[1].其缺點是需要完備的投影數(shù)據(jù).如稀疏視角下的FBP重建算法會引起條狀偽影[2].隨著壓縮感知等稀疏優(yōu)化理論的提出,基于優(yōu)化的迭代算法成為了十余年來的研究熱點[3,4].其中,總變差最小(total variation minimization,TV)算法是最成功的優(yōu)化算法之一,其已經在平板CT[5]、偏置掃描平板CT[6]、短掃描平板CT[7],C-arm CT[8]、低劑量平板CT[9]及機載成像儀(on-board imager,OBI)平板CT等[10]三維CT中得到了成功的應用,展示了其高精度的稀疏重建能力.

    TV算法可以高精度地稀疏重建的原因可以從多個方面來解釋.TV算法體現(xiàn)了壓縮感知的思想.壓縮感知認為,如果物體在某個變換域是稀疏的,那么可以僅僅使用稀疏采集到的原始數(shù)據(jù)來高精度地重建物體.其實現(xiàn)策略是在數(shù)據(jù)保真的約束下,使得物體的稀疏變換的?1范數(shù)最小[11].TV算法使用的稀疏變換是梯度大小變換(gradient magnitude transform,GMT)[12].TV算法體現(xiàn)了先驗約束的思想.在稀疏重建時,迭代法構建的線性方程組是欠定的,使得該線性方程組有無窮多個解.TV算法可以被認為加入了一種先驗知識,即物體的GMT是稀疏的.這樣,TV算法是從無窮多個解中選取TV最小的解.TV算法也體現(xiàn)了低通濾波的作用.TV算法最初被當作去噪方法應用到圖像處理中,其可以在去噪的同時保持圖像的邊緣信息[13].FBP的稀疏重建引入的條狀偽影是一種高頻噪聲.TV正是通過對這種噪聲的低通濾波,消除了條狀偽影,實現(xiàn)了高精度稀疏重建.

    TV最優(yōu)化模型有約束形式和非約束形式兩種類型.非約束形式的模型參數(shù)是一個平衡因子,用來平衡數(shù)據(jù)保真力度和正則力度,其沒有明確的物理意義,難以選擇.約束形式的TV模型,其模型參數(shù)有物理意義,更容易選擇[9].常用的約束的TV模型,以數(shù)據(jù)保真項為約束,TV正則項為目標函數(shù),其缺點是數(shù)據(jù)保真項對應的數(shù)據(jù)容差限這一參數(shù)對重建質量敏感,選取方法魯棒性較差[12].芝加哥大學Pan研究組從2014年陸續(xù)提出了TV約束的、數(shù)據(jù)分離最小(TV constrained,data divergence minimization,TVcDM)新型TV模型,并將之應用到了扇形束ROI重建[14],C-arm CT[8]、短掃描平板CT[7]及正電子發(fā)射斷層掃描(positron emission computed tomography,PET)[15]中.該模型的模型參數(shù)是TV限,其在較大范圍的變化對圖像質量影響較小,故易于選取,選取方法魯棒性較好.

    TVcDM模型是一種約束的最優(yōu)化模型,傳統(tǒng)的非約束優(yōu)化模型求解算法,如交替方向乘子法(alternating direction method of multipliers,ADMM)和分裂布萊格曼法(split Bregman,SB)[16],不適宜用來求解該模型.Chambolle-Pock(CP)算法框架[17?20]是一種優(yōu)秀的最優(yōu)化算法框架,不但適用于非約束的優(yōu)化模型,而且適用于約束的優(yōu)化模型.它有諸多優(yōu)點:其一,算法參數(shù)均可以計算得到而不需人為設定;其二,可以解所有凸最優(yōu)化模型,無論模型平滑還是非平滑;其三,算法的每個子最優(yōu)化問題一般均有閉合解.CP算法框架可以用來推導求解TVcDM模型的CP算法實例.

    Pan研究組[20]于2012年推導了一系列CP算法實例.這些實例分別用來求解最小二乘模型、非約束的?2-TV模型、非約束的?1-TV模型、非約束的KL-TV模型以及數(shù)據(jù)分離約束的、TV最小(data divergence constrained,TV minimization,DDcTV)模型,但是當時尚未提出TVcDM模型,所以并沒有給出TVcDM模型的CP算法實例的推導.2014—2017年,該組的一系列使用TVcDM模型的CT和PET圖像重建的文章,主要聚焦在完全重建、有限角度或者投影被截斷情況下的實際重建問題.2018年,本文作者與Pan研究組及芝加哥大學EPRI中心研究人員將這一模型應用于電子順磁共振成像(electron paramagnetic resonance imaging,EPRI)的稀疏重建中[21].這些研究聚焦在模型和算法在實際數(shù)據(jù)中的評估,而沒有使用仿真數(shù)據(jù),即存在金標準的情況下,從定量的角度展開收斂規(guī)律、稀疏重建能力評估和參數(shù)選擇等理論和方法層面的研究.

    鑒于如上原因,本文在總結出CP算法實例推導方法的基礎上,詳細推導了TVcDM-CP算法實例,并給出算法的詳細解釋;以Shepp-Logan和FORBILD仿真模體[22]的CT重建為例,實現(xiàn)該算法,驗證模型和算法的正確性;通過分析多個度量標準的迭代行為,揭示其收斂規(guī)律;通過稀疏重建,評估該模型的高精度稀疏重建能力;通過對含噪投影數(shù)據(jù)的重建,研究了模型參數(shù)對重建質量的影響.最后研究了算法參數(shù)對收斂速率的影響.

    2 方 法

    連續(xù)-連續(xù)(continuous-to-continuous,CC)成像模型把投影和圖像均看成連續(xù)函數(shù).而離散-離散(discrete-to-discrete,DD)成像模型把投影和圖像均看成離散函數(shù).這樣,DD模型就是一個線性方程組,基于DD模型的重建就是求方程組的解.一般情況下,該線性方程組是一個大規(guī)模的、病態(tài)的且往往欠定的(如稀疏重建時)方程組,使得求解成為需要重點攻克的問題.為了得到所需要的有用解,往往需要將此成像模型進一步建模為一個最優(yōu)化問題,使其可以將稀疏先驗等先驗知識引入模型.其后,需要設計求解算法,以解最優(yōu)化模型.

    基于如上脈絡,以CT重建為背景,以TVcDM模型為研究對象,在該部分將按照成像系統(tǒng)建模-最優(yōu)化模型-CP算法-參數(shù)選擇這一研究鏈展開敘述.因為本文旨在探索新穎的TV模型及其求解算法,不失一般性,以平行束CT作為研究對象.別的成像模態(tài)或者CT成像的其他掃描模式遵循同樣的規(guī)律,其研究方法只需要做相應的類推或擴展.

    2.1 成像模型

    平行束CT的DD成像模型是一個線性方程組,如(1)式所示:

    這里,g是一個大小為M的列向量,表示離散的投影數(shù)據(jù);u是一個大小為N的列向量,表示離散圖像;A是一個大小為M×N的矩陣,這里代表2D Radon變換;Am,n表示第n個像素對第m個投影測量值的貢獻.對2D Radon變換來說,Am,n是第n個像素(正方形)與第m條射線(直線)的相交長度.

    設圖像的大小為nx×ny,將其一列一列串行化,就可以得到圖像的向量表示,其大小為N=nx×ny.設投影集大小為np×na,即共有na個角度下的投影,每個投影上有np個點,將其一列一列串行化,就可以得到投影數(shù)據(jù)的向量表示,其大小為M=np×na.

    實際上,求取Am,n的方法就是所謂的投影方法,因為投影是系統(tǒng)矩陣與圖像向量的乘積.投影方法包括像素驅動法[23],Siddon射線驅動法[24],Joseph射線驅動法[25]和距離驅動法[26]等.傳統(tǒng)的像素驅動法會引入高頻震蕩偽影,已經設計了一種精確的像素驅動法[23];本文使用了精確的像素驅動法來生成投影,即求取系統(tǒng)矩陣.

    2.2 最優(yōu)化模型

    由于(1)式所示的線性方程組的病態(tài)、大規(guī)模及欠定等因素,需要進一步將之建模為一個最優(yōu)化模型.根據(jù)重建的需要,可以將壓縮感知、低秩矩陣等理論及各種先驗知識耦合,以設計最優(yōu)化模型.本文研究一種新穎的基于壓縮感知的TV模型——TVcDM.

    TVcDM,即TV約束的、數(shù)據(jù)分離最小模型,可以表示為

    這里,u?是最優(yōu)化模型的解,即被重建的圖像;∥.∥表示?2范數(shù)的平方;圖像的TV范數(shù)∥u∥TV是圖像梯度大小變換|Du|mag的?1范數(shù):

    這里D是一個大小為2N ×N的矩陣

    D1和D2均為N×N的矩陣,分別表示沿著x和y軸方向的離散梯度變換.

    令ux,y,x∈[1,nx],y∈[1,ny]表示二維圖像的每個像素,則D1變換可以表示為

    D2變換可以表示為

    令|.|mag表示對一個二維向量求模,可以是?1范數(shù)模,也可以是?2范數(shù)模.在本文中,采用?2范數(shù)模,即

    這樣,(3)式的各向同性形式可以表示為

    為了提高收斂速度,引入λ和ν兩個算法參數(shù)以調整TV范數(shù)和數(shù)據(jù)保真項這兩個凸函數(shù)的相對大小,如此得到的新的TVcDM模型為[7,8,14]

    這兩個參數(shù)雖然出現(xiàn)在最優(yōu)化模型中,但是這里稱其為算法參數(shù),因為它們不能決定最優(yōu)化模型的解,但是可以影響解的收斂路徑和速度.

    2.3 Chambolle-Pock算法及其解釋

    表1給出了用于求解最優(yōu)化模型(9)的CP算法偽代碼.

    在表1中,∥.∥sv是矩陣的最大奇異值,其求法見文獻[20]中的算法8;向量u,uˉ及s的大小為N;向量p的大小為M;向量q和a的大小為2N;6.3中的向量1I是一個大小為N 的“1”向量;6.2中的投影操作符ProjectOnto?1Ballνt1的求法見表2的偽代碼;uconv是達到設定的實用收斂條件后的收斂解.

    在表2中,x表示一個長度為N的向量;a是?1球的半徑;m是一個長度為N的向量,其每個元素對應x中相應元素的絕對值;sign(x)是符號函數(shù):正數(shù)的函數(shù)值是1,0的函數(shù)值是0,負數(shù)的函數(shù)值是?1.

    表1 用于求解最優(yōu)化模型(9)的CP算法偽代碼[7,8,14]Table 1.Pseudocode of the CP algorithm for solving the optimization model shown in(9)[7,8,14].

    表2 ProjectOnto?1Ball的求解算法偽代碼[14]Table 2. Pseudocode of the solving algorithm for ProjectOnto?1Ball[14].

    2.4 CP算法的推導步驟

    本節(jié)首先從CP算法框架總結出算法實例的推導方法,然后推導TVcDM模型的CP算法實例.

    2.4.1 CP算法框架及其推導方法

    CP算法框架可以解形如(10)式所示的最優(yōu)化模型.

    在此模型中,x和y是分別屬于向量空間X和Y的多維向量;K是一個矩陣,表示一個線性變換,定義了一個從X到Y的線性映射;G和F分別是關于x和y的多元函數(shù),分別定義了從X和Y到實數(shù)的一個映射.CP框架要求這兩個函數(shù)是凸函數(shù),但不要求平滑.

    CP算法框架見表3.

    表3 CP算法框架(N步迭代)[17]Table 3.The CP algorithm framework(N steps of iterations)[17].

    表3中,∥K∥SV是指矩陣K的最大奇異值;F?是指凸函數(shù)F 的共軛凸函數(shù);proxσ[F?]和proxτ[G]是最鄰近映射操作;T表示轉置.

    任意一個凸函數(shù)H(z)的共軛凸函數(shù)定義為

    這里,需要注意,max操作符不是求取使得目標函數(shù)最大的自變量,而是求取當c為最大化值時目標函數(shù)的值.

    prox最鄰近映射的表達式為

    該最鄰近映射本質上是一個最優(yōu)化問題,但是因為z是一個變量,所以其結果是一個關于z的函數(shù).

    現(xiàn)在,可以總結出使用CP算法框架推導CP算法實例的步驟為:

    1)構造如(10)式所示的最優(yōu)化問題;

    2)求取F?;

    3)求 取proxσ[F?]和proxτ[G]兩 個 最 鄰 近映射;

    4)代入表3的第4行和第5行形成算法實例.

    在如上步驟中,關鍵問題是基于(11)式的凸共軛函數(shù)的求取和基于(12)式的最鄰近映射的求取.CP算法的優(yōu)勢之一是每個子最優(yōu)化問題都有閉合解.這就要求在(12)式中,函數(shù)H(Z)要足夠簡單,以使得可以求得閉合形式的解.現(xiàn)有研究表明,在圖像重建中使用的大部分凸函數(shù)所對應的最鄰近映射均有閉合解[20].

    2.4.2 TVcDM-CP算法實例的推導

    (9)式所對應的最優(yōu)化問題,可以通過如下變換形成(10)式的形式.

    其中,δ?1Ball是一個指示函數(shù),表示如果一個n維向量在特定半徑的?1球內部,則函數(shù)值為0,否則,值為∞,其表達式為:該式表示,如果向量x在半徑為a的?1球內部,則函數(shù)值為0,否則值為∞.

    現(xiàn)在,根據(jù)如下所示的一系列變換,將(12)與(10)式對應起來.

    這樣,

    根據(jù)CP算法實例的推導步驟,現(xiàn)在的任務是求取兩個凸函數(shù)F1(p)和F2(q)的共軛凸函數(shù),以及求取FF和G三個函數(shù)的prox最鄰近映射.

    根據(jù)(11)式,可得

    根據(jù)(12)式,可得

    將(15),(19)—(21)代入表3,可以得到表1所示的TVcDM模型的CP算法實例.

    2.5 重建參數(shù)

    最優(yōu)化模型及其算法構成了一套完整的重建方案,其中所有的參數(shù)構成了重建參數(shù).與模型相關的參數(shù)稱為模型參數(shù);與算法相關的參數(shù)稱為算法參數(shù).

    模型參數(shù)包括TV限t1、像素大小、投影方法等.在本文中,像素大小和探元大小都歸一化為1,投影方法采用精確的像素驅動法,故將在3.3節(jié)重點研究TV限對重建的影響.

    算法參數(shù)包括σ,τ,λ,ν和θ,因為σ,τ和θ參數(shù)均可以計算得到,而不需要人為設定,所以不需要關注.λ和ν對算法的收斂速度和路徑有影響,故將在3.4節(jié)研究其對收斂速率的影響.

    3 結 果

    在本節(jié)中,設計4個研究點:1)重建方法的驗證及CP算法的收斂行為分析;2)重建方法的稀疏重建能力評估;3)TV限對重建精度的影響;4)λ和ν的選擇對收斂速率的影響.

    3.1 逆犯罪研究——重建方法的驗證及CP算法的收斂行為

    從數(shù)學的角度講,圖像重建是一個反問題.對于平行束CT重建,解析法,如濾波反投影算法,是反Radon變換問題;迭代法是線性方程組的求解問題.逆犯罪是一種評估反問題是否正確的驗證方法.當仿真數(shù)據(jù)理想且充分的情況下,如果求得的解的誤差在計算機精度范圍內可以充分小,則認為逆犯罪成功,表明設計的整套重建方案及其實現(xiàn)——從成像系統(tǒng)建模、最優(yōu)化問題建模、算法設計到計算機實現(xiàn),都是正確的.

    在本部分中,采用Shepp-Logan和FORBILD模型為仿真模型;仿真生成理想而充分的投影數(shù)據(jù);使用第二部分的重建方法進行重建;通過逆犯罪研究驗證整套重建方法的正確性.

    3.1.1 仿真模體及重建條件

    分別采用Shepp-Logan模體和FORBILD模體為仿真模型,圖像大小為256×256,每個像素被看成單位1大小的像素.旋轉中心定位于128×128.探測器長度為256,每個探元大小為單位1.均勻采集[0,π]范圍360個角度下的投影,所以投影數(shù)據(jù)的大小為256×360.采用精確的像素驅動法定義系統(tǒng)矩陣A,并生成投影數(shù)據(jù).如此構建的線性方程組,有256×256個未知數(shù),有256×360個方程,屬于過定方程.但是因為方程組是相容的,所以方程組有惟一的精確解.

    在CP算法中,設定λ=1,b=1即ν=νA,其余算法參數(shù)按照表1中計算得到.TV限設定為仿真模型的TV值.

    CP迭代算法的收斂條件定義為RMSE(un,utruth)6 10?4.RMSE的定義見(22)和(23)式.

    3.1.2 重建結果的定性及定量分析

    重建結果采用定性觀察和定量分析的方法.一方面視覺檢查圖像質量,一方面比較圖像中線位置的波形,以進行定性觀察.定量分析,則采用均方根誤差(root-mean-square-error,RMSE)作為度量標準,如下式所示:

    式中x是待評估向量,s是向量真值,n是向量的長度.

    對Shepp-Logan重建來說,CP算法在迭代到2273次時,達到了設定的收斂條件.重建結果見圖1. 圖1(a)是Shepp-Logan仿真模型圖像;圖1(b)是收斂態(tài)的重建圖像.比較這兩個圖像可以發(fā)現(xiàn)它們幾乎完全一樣,用肉眼難以分辨彼此.圖1(c)是重建圖像和仿真模型圖像的中心線波形的比較,從圖中可以看出,兩條波形幾乎完全重合.圖像和波形的主觀評價表明取得了高精度重建.

    對FORBILD重建來說,CP算法在迭代到1712次時,達到了設定的收斂條件.重建結果見圖2.圖2(a)是FORBILD仿真模型圖像;圖2(b)是收斂態(tài)的重建圖像.比較這兩個圖像可以發(fā)現(xiàn),它們幾乎完全一樣,用肉眼難以分辨彼此;圖2(c)是重建圖像和仿真模型圖像的中心線波形比較圖.從圖中可以看出,兩條波形幾乎完全重合.圖像和波形的主觀評價表明取得了高精度重建.

    兩種模體的重建圖像的RMSE均達到了10?4,一方面表示算法的收斂條件達到了,另一方面表明逆犯罪成功,因為此反問題被高精度地求解了出來.

    圖1 Shepp-Logan仿真模體的TVcDM重建結果 (a)Shepp-Logan仿真模體圖像;(b)重建圖像;(c)模型圖像和重建圖像左右中心線波形的比較:紅色線是真實波形,藍色線是重建波形,二者幾乎完全重合;圖像的顯示窗口是[0,1]Fig.1.The TVcDM reconstruction images of the Shepp-Logan simulation phantom:(a)The Shepp-Logan phantom image;(b)the reconstructed image;(c)comparison of the profiles on the vertical center-line of the phantom image and the reconstructed image.In(c),it may be seen that the red truth profile and the blue reconstructed profile are almost completely overlapped.

    圖2 FORBILD仿真模體的TVcDM重建結果 (a)FORBILD仿真模體圖像;(b)重建圖像;(c)模體圖像和重建圖像左右中心線波形的比較圖:紅色線是真實波形,藍色線是重建波形,二者幾乎完全重合;圖像的顯示窗口是[0,1]Fig.2.The TVcDM reconstruction images of the fORBILD simulation phantom:(a)the fORBILD phantom image;(b)the reconstructed image;(c)comparison of the profiles on the vertical center-line of the phantom image and the reconstructed image.In(c),it may be seen that the red truth profile and the blue reconstructed profile are almost completely overlapped.

    3.1.3 CP算法的收斂行為

    為了描述收斂行為,引入3個度量標準,通過觀察3個度量標準的迭代走勢,揭示CP算法的迭代規(guī)律.這些度量標準如下式所示:

    (23)式描述的是重建的圖像與仿真模型圖像的RMSE距離;(24)式描述的是重建圖像生成的投影與原始投影之間的?2范數(shù)距離;(25)式描述重建圖像的TV值與TV真值的相對誤差.在此逆犯罪研究中,因為模型的投影和圖像是完全一致的、數(shù)據(jù)是充分且精確的,所以這些度量標準在收斂態(tài)均可以充分小,表示取得了高精度重建.

    圖3和圖4分別展示了Shepp-Logan和FORBILD模體的CP算法的收斂行為,圖3(a)—(c)分別展示了(23)—(25)式所示的度量標準的迭代走勢.

    從圖3和圖4的(a)可以看出,圖像誤差不斷下降,達到了設定的收斂條件,而且可以看到其有繼續(xù)下降的趨勢.如果繼續(xù)推進迭代,圖像誤差甚至可以達到10?6.灰度圖像的灰度只有256個等級,即從0到255,理論上說,當圖像誤差達到1/256≈3.9‰時,顯示器已經無法分辨存在此誤差水平的兩幅圖像.所以,我們設定逆犯罪標志,即收斂條件為圖像誤差小于等于10?4.

    圖3和圖4的(b)顯示的是數(shù)據(jù)誤差,即數(shù)據(jù)殘差的走勢.數(shù)據(jù)殘差是TVcDM模型的目標函數(shù).可以看出,數(shù)據(jù)殘差不斷下降,而且在當前的停止點有繼續(xù)向下的趨勢.

    圖3 TVcDM-CP算法重建Shepp-Logan模體時的收斂行為 (a)重建誤差M1(n)的迭代走勢;(b)數(shù)據(jù)誤差M2(n)的迭代走勢;(c)重建圖像的TV相對誤差M3(n)的迭代走勢.Fig.3.The convergence behavior of the TVcDM-CP algorithm for reconstructing the Shepp-Logan phantom:(a)The iteration trend of the reconstruction error M1(n);(b)the iteration trend of the data error M2(n);(c)the iteration trend of the relative TV error M3(n).

    圖4 TVcDM-CP算法重建FORBILD模體時的收斂行為 (a)重建誤差M1(n)的迭代走勢;(b)數(shù)據(jù)誤差M2(n)的迭代走勢;(c)重建圖像的TV相對誤差M3(n)的迭代走勢Fig.4.The convergence behavior of the TVcDM-CP algorithm for reconstructing the fORBILD phantom:(a)The iteration trend of the reconstruction error M1(n);(b)the iteration trend of the data error M2(n);(c)the iteration trend of the relative TV error M3(n).

    圖3和圖4的(c)顯示了TV值相對于TV限的相對誤差,其整體趨勢向下,但是在兩個重建案例中,當?shù)胶笃?出現(xiàn)了較大抖動,相對誤差變大.這種現(xiàn)象是該算法迭代行為的一個特點,在我們各種各樣的重建實踐中經常出現(xiàn),這并不影響算法的收斂性.雖然TV值相對誤差出現(xiàn)了向上抖動,但是圖像誤差和數(shù)據(jù)誤差是整體不斷下降的,表明最優(yōu)化模型的解仍然在向收斂態(tài)逼近.其實,可以觀察到TV值下降到了一個很小的值,比如10?7之后產生了向上的抖動,其達到一定值之后,又會繼續(xù)向下.此后,即使再有抖動,其上限也會保持在一個很小的值附近,比如10?4.而TV值是指圖像梯度變換的?1范數(shù),TV值在此極小范圍的波動,分配到圖像每個像素上(在本節(jié)案例中,圖像包含256×256=65536個像素),圖像基本沒有變化.

    從圖3和圖4的所有子圖可以看出,所有的迭代曲線均存在抖動現(xiàn)象,這是TVcDM-CP算法迭代行為的一個特點.圖像誤差和數(shù)據(jù)殘差(目標函數(shù))整體會不斷下降,直到計算機浮點數(shù)精度(本文因為投影與反投影操作使用GPU加速,所以均采用了單精度數(shù))所允許的收斂態(tài).

    3.2 重建方法的稀疏重建能力評估

    為了評估TVcDM模型的稀疏重建能力,在仿真實驗中,分別從30,60,90和120個角度重建圖像,分析投影個數(shù)對重建的影響.重建條件與3.1節(jié)相同,除了Shepp-Logan模體重建的收斂條件變?yōu)榈?000次結束,而FORBILD模體重建的收斂條件為迭代至2000次結束.

    圖5是Shepp-Logan模體的TVcDM的重建結果與FBP算法重建結果的比較.由圖可以看出隨著投影個數(shù)從120個減少到30個,FBP重建的條狀偽影越來越嚴重,30個角度下的FBP重建引入了明顯的偽影;同時看到,各個角度下的TVcDM重建圖像幾乎完全相同,看不出明顯的因為投影個數(shù)減少而產生的質量退化現(xiàn)象.30個角度下的投影已經可以實現(xiàn)TVcDM的高精度重建,表明TVcDM模型有高精度稀疏重建能力.

    圖5 Shepp-Logan模體重建中TVcDM與FBP算法的重建結果比較 圖像上面的數(shù)字表示投影個數(shù),左邊的文字表示所使用的算法;顯示窗口為[0,1]Fig.5.Comparison of the reconstructed Shepp-Logan images by use of TVcDM and FBP algorithm.The numbers above the image indicate projection number and the text at the left of the image indicate the algorithm used.The display window is[0,1].

    圖6是FORBILD模體的TVcDM的重建結果與FBP算法重建結果的比較.由圖可以看出,隨著投影個數(shù)從120個減少到30個,FBP重建圖像的條狀偽影越來越嚴重,60個角度下的FBP重建引入了明顯的偽影而30個角度下的偽影已經異常嚴重;同時看到,各個角度下的TVcDM重建圖像幾乎完全相同,看不出明顯的因為投影個數(shù)減少而產生的質量退化現(xiàn)象.30個角度下的投影已經可以實現(xiàn)TVcDM的高精度重建,表明TVcDM模型有高精度稀疏重建能力.

    為了定量比較TVcDM與FBP算法在不同投影個數(shù)情形下的重建精度,圖7繪制了重建圖像的RMSE相對投影個數(shù)的走勢圖.圖7(a)對應Shepp-Logan模體重建;圖7(b)對應FORBILD模體重建.兩個模體重建對應的RMSE走勢規(guī)律是相同的.隨著投影個數(shù)從120減小到30,FBP重建結果的RMSE逐漸增大,表示FBP算法對投影個數(shù)的變化是敏感的,其沒有高精度稀疏重建能力,當投影個數(shù)為30個時,其RMSE過大,圖像精度過低,導致重建圖像失去有效性.而TVcDM算法對投影個數(shù)變化不敏感,投影個數(shù)為60,90和120時,RMSE非常接近.投影個數(shù)為30時,因為投影視角過于稀疏,重建誤差增大.但即使TVcDM算法只使用30個角度下的投影,其重建圖像的精度也遠遠高于FBP算法從120個角度下的投影重建圖像的精度.

    圖6 FORBILD模體重建中TVcDM與FBP算法的重建結果比較 圖像上面的數(shù)字表示投影個數(shù),左邊的文字表示所使用的算法;顯示窗口為[0,1]Fig.6.Comparison of the reconstructed FORBILD images by use of TVcDM and FBP algorithm.The numbers above the image indicate projection number and the text at the left of the image indicate the algorithm used.The display window is[0,1].

    圖7 重建圖像RMSE隨投影個數(shù)變化的走勢圖 (a)對應Shepp-Logan模體重建;(b)對應FORBILD模體重建Fig.7.Plot of RMSE of the reconstructed image as function of projection number with(a)corresponding to Shepp-Logan phantom reconstruction and(b)FORBILD phantom reconstruction.

    圖5和圖6的定性主觀視覺檢查和圖7的定量誤差分析均表明:TVcDM-CP算法有從稀疏投影中高精度重建圖像的能力.圖7(b)中30個角度對應的重建RMSE比圖7(a)中30個角度對應的重建RMSE大.此現(xiàn)象表明:不同的重建對象,因為復雜程度不同,所需的最少投影個數(shù)是不同的.FORBILD圖像比Shepp-Logan圖像復雜,所以要達到特定的重建精度,需要的投影個數(shù)更多.

    3.3 TV限對重建精度的影響

    在仿真實驗中,給投影加45.0 dB的高斯白噪聲,進行重建.投影個數(shù)為30個;TV限分別設定為TV真值的0.6,0.8,1.0,1.2和1.4倍;其余重建條件相同.

    圖8和圖9是Shepp-Logan模體重建結果.圖8是不同TV限約束下的對含噪投影的重建圖像.圖9是不同TV限約束下的對含噪投影的重建圖像與仿真模體(真值)的差分絕對值圖像.從圖8(c)和圖9(c)可以看出,當TV限選取正確,即選取TV真值時,可以獲得高精度的重建結果,TVcDM模型不但壓制了因為稀疏采集(只有30個投影)可能引起的條狀偽影,而且壓制了因為高斯白噪聲污染而產生的圖像噪聲.從圖8和圖9的(a)和(b)可以看出,當TV限選取過小時,TV范數(shù)的平滑作用過大,使得重建圖像產生了低頻過平滑偽影.而從圖8和圖9的(d)和(e)可以看出,當TV限選取過大時,TV范數(shù)的平滑作用減小,使得重建圖像產生了高頻噪聲.可見,TV限的選擇對重建有重要影響,其決定了最優(yōu)化模型的收斂解.只有在平滑和細節(jié)保持之間,尋找到平衡點,選取出最適合的TV限,才可以高精度地重建圖像.

    圖10(a)是Shepp-Logan模體重建中RMSE隨TV限變化的走勢圖.可以看到TV限過大或者過小,都導致了更大的誤差,選取合適的TV限會取得較小誤差,實現(xiàn)高精度重建.

    圖11和圖12是FORBILD模體重建結果.結合圖10(b),可以看出,過小或者過大的TV限均導致了更大的誤差,合適的TV限可以取得較小誤差,取得高精度重建.

    圖8 TV限的不同選擇對重建精度的影響 (a)—(e)分別是TV限為TV真值的0.6,0.8,1.0,1.2和1.4倍時對應的重建圖像;為了突出顯示,采用了偽彩色顯示方式;顯示窗口是[0,1]Fig.8.The impact of the selection of TV tolerance on reconstruction accuracy:(a)–(e)are the reconstructed images with TV tolerance being 0.6,0.8,1.0,1.2 and 1.4 times of the TV real value.To highlight the display of the images,the images are displayed in pseudocolor pattern and the display window is[0,1].

    圖9 TV限的不同選擇對重建精度的影響 (a)—(e)分別是TV限為TV真值的0.6,0.8,1.0,1.2和1.4倍時對應的重建圖像與仿真模體的差分絕對值圖像;為了突出顯示,采用了偽彩色顯示方式;顯示窗口是[0,0.2]Fig.9.The impact of the selection of TV tolerance on reconstruction accuracy:(a)–(e)are the reconstructed absolute-difference images with TV tolerance being 0.6,0.8,1.0,1.2 and 1.4 times of the TV real value.To highlight the display of the images,the images are displayed in pseudocolor pattern and the display window is[0,0.2].

    兩種模體復雜程度不同,圖像特征不同,但是在此研究中,TV限對重建精度的影響規(guī)律一致.實驗結果和定性分析一致表明,TV限對重建精度有重要影響,只有選擇最佳TV限才可以在圖像的平滑度和細節(jié)保持之間找到最佳平衡點,實現(xiàn)高精度重建.

    圖10 RMSE隨TV限變化的走勢圖 (a)和(b)分別對應Shepp-Logan重建和FORBILD重建Fig.10.Plot of RMSE as function of TV tolerance with(a)corresponding to Shepp-Logan reconstruction and(b)FORBILD reconstruction.

    圖11 TV限的不同選擇對重建精度的影響 (a)—(e)分別是TV限為TV真值的0.6,0.8,1.0,1.2和1.4倍時對應的重建圖像;為了突出顯示,采用了偽彩色顯示方式;顯示窗口是[0,1]Fig.11.The impact of the selection of TV tolerance on reconstruction accuracy:(a)–(e)are the reconstructed images with TV tolerance being 0.6,0.8,1.0,1.2 and 1.4 times of the TV real value.To highlight the display of the images,the images are displayed in pseudocolor pattern and the display window is[0,1].

    圖12 TV限的不同選擇對重建精度的影響 (a)—(e)分別是TV限為TV真值的0.6,0.8,1.0,1.2和1.4倍時對應的重建圖像與仿真模體的差分絕對值圖像;為了突出顯示,采用了偽彩色顯示方式;顯示窗口是[0,0.2]Fig.12.The impact of the selection of TV tolerance on reconstruction accuracy.(a)–(e)are the reconstructed absolute-difference images with TV tolerance being 0.6,0.8,1.0,1.2 and 1.4 times of the TV real value.To highlight the display of the images,the images are displayed in pseudocolor pattern and the display window is[0,0.2].

    3.4 λ和ν的選擇對收斂速率的影響

    仍然使用3.1.1中設定的重建條件,只是將投影個數(shù)變?yōu)?0.采用Shepp-Logan和FORBILD兩個模體分別進行實驗.為了分析λ和ν的不同選擇對重建的影響,設計5種不同的組合形式,使算法均迭代至1000次,觀察RMSE的下降態(tài)勢,比較其對重建速率的影響.

    圖13是收斂速率比較圖,其中圖13(a)對應Shepp-Logan重建;圖13(b)對應FORBILD重建.

    從圖13可以看出,不同的λ和ν對應不同的收斂速率.對于兩個模體重建來說,λ=1,ν=0.1νA均為收斂速率最快的一個組合.但是如果按照從慢到快排序,對兩個模體重建而言,順序不完全相同.λ和ν的最佳選擇是依賴于特定重建對象的.即使對于同一個重建對象,當投影個數(shù)發(fā)生變化時,λ和ν的最快組合也可能發(fā)生變化.對于一個特定的重建任務,要想找到最快的組合,需要使用不同的組合重建,通過觀察迭代收斂行為,選擇最快的一個組合.

    圖13 不同的λ和ν的組合下RMSE在1000次迭代范圍內的收斂速率比較 (a)和(b)分別對應Shepp-Logan重建和FORBILD重建Fig.13. Comparison of the convergence rates of RMSE of reconstructions with different λ and ν in the iteration range of 1 to 1000 with(a)and(b)corresponding to Shepp-Logan and FORBILD reconstruction,respectively.

    4 討 論

    4.1 模型設計相關問題

    為了提高TV模型的精確重建能力,人們提出了各種TV范數(shù)的變種:廣義的各向異性TV范數(shù),保邊TV范數(shù),適用于有限角度重建的各向異性TV范數(shù),高階TV范數(shù);TpV(p∈[0,1])范數(shù)等.類似的最優(yōu)化模型還包括小波變換最小、基于字典學習的稀疏模型以及非局部均值等.所有這些模型均可以用一個統(tǒng)一優(yōu)化模型來描述.在統(tǒng)一優(yōu)化模型具體化時,其一要考慮選取何種正則模型,其二要考慮選取何種數(shù)據(jù)保真模型,其三則是考慮選擇約束形式還是非約束形式.

    最優(yōu)化模型是數(shù)學模型,它是對物理模型的一種近似,因此各種各樣的最優(yōu)化模型沒有優(yōu)劣之分,只有在特定重建場合下的適用性問題.對于圖像重建的一個特定任務,應用驅動地設計最適合的最優(yōu)化模型是基于優(yōu)化的重建方案設計的最重要的工作.

    最優(yōu)化模型決定了最終的解,所以模型參數(shù)對重建有重要影響.當前,盡管人們已經系統(tǒng)研究了模型參數(shù)對重建的影響,但是模型參數(shù)的自動選擇問題應該是今后的研究重點.

    4.2 關于重建算法

    本文所用模型的求解是非常具有挑戰(zhàn)性的:1)規(guī)模巨大,系統(tǒng)矩陣很大,設計的算法中不能出現(xiàn)矩陣求逆操作;2)往往欠定,方程個數(shù)遠小于未知量個數(shù);3)病態(tài),因為噪聲和數(shù)據(jù)的不一致性,造成方程組是不相容的;4)TV范數(shù)是一個不平滑的凸函數(shù).

    人們已經為各種各樣的最優(yōu)化模型設計了求解算法,如ASD-POCS算法[12],ADMM算法以及本文使用的CP算法等.需要注意的是,模型決定解,求解算法僅僅決定解的收斂路徑和速度.所以,各種各樣的求解算法的重建精度是沒有可比性的,惟一可以比較的是其收斂速度.本文采用CP算法的原因是該算法保證收斂、算法參數(shù)不需要調節(jié)且每個子最優(yōu)化問題均有閉合解.

    5 結 論

    本文設計了TVcDM-CP重建方法.將圖像重建問題建模為一個TVcDM最優(yōu)化模型.在該約束的最優(yōu)化模型中,TV正則項是約束項,數(shù)據(jù)保真項是目標函數(shù).基于CP算法框架,詳細推導了TVcDM模型的CP算法實例;經仿真模型重建,驗證了模型及算法的正確性并分析了算法的收斂行為;評估了模型的稀疏重建能力;分析了模型參數(shù)的選擇對重建的影響.

    TVcDM模型相較于DDcTV模型,一個可能的優(yōu)勢是模型參數(shù)比較容易選擇.DDcTV模型的模型參數(shù)是數(shù)據(jù)容差限[12],其大小決定了重建圖像的投影和實際的測量投影的距離.在實際重建中,投影數(shù)據(jù)包含噪聲和不一致性,只有選擇合適的數(shù)據(jù)容差限才可以有效抑制噪聲和不一致性.選取過小,則噪聲被引入;選取過大則重建圖像過于平滑,細節(jié)丟失.TVcDM的模型參數(shù)是TV限,合理的TV限估計會實現(xiàn)高精度重建,而選取過小則圖像過于平滑,選取過大則圖像被引入噪聲.TV限在小范圍內的變化對重建影響不大(參考文獻[8]的圖7),所以該參數(shù)的選擇比較魯棒.

    本文研究表明,TVcDM-CP重建方法可以高精度地從稀疏角度下的投影中重建圖像;在CP算法的收斂過程中,會出現(xiàn)抖動現(xiàn)象,但是整體呈現(xiàn)下降趨勢;TV限對重建質量有重要影響,合理選取該參數(shù)可以取得高精度的重建質量;λ和ν的不同取值組合會導致不同的收斂速率.

    本文所設計的TVcDM-CP算法是一種通用的圖像重建模型,可以被應用到平行束CT、扇形束CT、錐形束CT中,也可以被應用到PET,MRI及EPRI中.依據(jù)本文的研究鏈,即成像系統(tǒng)建模-最優(yōu)化問題建模-CP算法設計-算法實現(xiàn)及驗證-實際重建,各種重建模態(tài)均可以實現(xiàn)高精度稀疏重建.在此研究鏈中,需要做的改變僅僅是將本文的成像模型中的系統(tǒng)矩陣變?yōu)橄鄳某上衲B(tài)對應的系統(tǒng)矩陣,然后任務驅動地選擇模型參數(shù)和算法參數(shù).

    猜你喜歡
    模體范數(shù)高精度
    基于Matrix Profile的時間序列變長模體挖掘
    植入(l, d)模體發(fā)現(xiàn)若干算法的實現(xiàn)與比較
    高抗擾高精度無人機著艦縱向飛行控制
    基于加權核范數(shù)與范數(shù)的魯棒主成分分析
    基于網(wǎng)絡模體特征攻擊的網(wǎng)絡抗毀性研究
    矩陣酉不變范數(shù)H?lder不等式及其應用
    船載高精度星敏感器安裝角的標定
    基于高精度測角的多面陣航測相機幾何拼接
    基于模體演化的時序鏈路預測方法
    自動化學報(2016年5期)2016-04-16 03:38:40
    高精度免熱處理45鋼的開發(fā)
    山東冶金(2015年5期)2015-12-10 03:27:41
    一区在线观看完整版| 国产精品免费大片| 亚洲色图综合在线观看| 岛国在线观看网站| videos熟女内射| 久久av网站| 亚洲av欧美aⅴ国产| 午夜老司机福利片| 久久青草综合色| 天堂俺去俺来也www色官网| 国产精品香港三级国产av潘金莲| 免费女性裸体啪啪无遮挡网站| 欧美日韩视频精品一区| 黑人猛操日本美女一级片| 午夜福利影视在线免费观看| 自线自在国产av| 国产欧美亚洲国产| 香蕉国产在线看| 制服人妻中文乱码| 成年人午夜在线观看视频| 亚洲成人免费av在线播放| 久久久国产精品麻豆| 久9热在线精品视频| 亚洲国产精品一区三区| 一级a爱视频在线免费观看| 亚洲精品一区蜜桃| 久久久国产精品麻豆| 国产欧美日韩综合在线一区二区| 久久久久久久精品精品| 动漫黄色视频在线观看| 欧美激情 高清一区二区三区| 国产精品 国内视频| 欧美一级毛片孕妇| 咕卡用的链子| 午夜免费成人在线视频| 91成年电影在线观看| 日韩有码中文字幕| 一二三四在线观看免费中文在| 老熟女久久久| 国产成人精品在线电影| 999久久久精品免费观看国产| 成年动漫av网址| 欧美xxⅹ黑人| 国产精品亚洲av一区麻豆| 欧美日韩国产mv在线观看视频| 日韩制服丝袜自拍偷拍| 久久精品人人爽人人爽视色| 日本一区二区免费在线视频| 脱女人内裤的视频| 国产精品一二三区在线看| 999精品在线视频| 久久精品国产亚洲av香蕉五月 | 亚洲欧美一区二区三区黑人| 99热国产这里只有精品6| 国产麻豆69| 女人精品久久久久毛片| 久久ye,这里只有精品| 精品卡一卡二卡四卡免费| 老汉色∧v一级毛片| 老司机靠b影院| 国产一区有黄有色的免费视频| 另类亚洲欧美激情| 国产精品.久久久| a级毛片在线看网站| 99久久精品国产亚洲精品| 两个人看的免费小视频| 啦啦啦视频在线资源免费观看| 男男h啪啪无遮挡| 国产成人免费无遮挡视频| 国产1区2区3区精品| 亚洲va日本ⅴa欧美va伊人久久 | 99国产综合亚洲精品| 免费不卡黄色视频| 9191精品国产免费久久| www.自偷自拍.com| 黑丝袜美女国产一区| 国产熟女午夜一区二区三区| 高潮久久久久久久久久久不卡| 欧美亚洲日本最大视频资源| 十八禁网站免费在线| 免费高清在线观看视频在线观看| 高清在线国产一区| 黑人巨大精品欧美一区二区mp4| 蜜桃在线观看..| 桃红色精品国产亚洲av| 亚洲精品国产色婷婷电影| 国产视频一区二区在线看| 欧美激情 高清一区二区三区| 嫁个100分男人电影在线观看| www.熟女人妻精品国产| 亚洲精品成人av观看孕妇| 久久中文看片网| 国产欧美日韩一区二区精品| 脱女人内裤的视频| 男人舔女人的私密视频| 国产精品自产拍在线观看55亚洲 | 超色免费av| 青春草亚洲视频在线观看| 亚洲欧美日韩另类电影网站| 少妇被粗大的猛进出69影院| 欧美精品一区二区大全| 亚洲欧洲日产国产| 精品福利观看| 亚洲精品久久成人aⅴ小说| 亚洲七黄色美女视频| 在线av久久热| 日韩免费高清中文字幕av| 国产野战对白在线观看| 中文字幕另类日韩欧美亚洲嫩草| 精品久久蜜臀av无| 欧美 亚洲 国产 日韩一| 国产亚洲av片在线观看秒播厂| 大陆偷拍与自拍| 国产区一区二久久| 大片免费播放器 马上看| 亚洲精品一卡2卡三卡4卡5卡 | 十八禁人妻一区二区| 亚洲精品乱久久久久久| 亚洲欧美一区二区三区久久| 极品少妇高潮喷水抽搐| 国产一区二区三区在线臀色熟女 | 中文欧美无线码| 国产色视频综合| 欧美日韩一级在线毛片| 大码成人一级视频| 黄色毛片三级朝国网站| 黄色视频在线播放观看不卡| 亚洲国产成人一精品久久久| 一区二区日韩欧美中文字幕| 人妻一区二区av| 亚洲精品中文字幕一二三四区 | 亚洲精品乱久久久久久| 久久久久久久精品精品| 欧美精品啪啪一区二区三区 | 99精品欧美一区二区三区四区| 在线观看免费日韩欧美大片| 国产成人av激情在线播放| 老熟妇仑乱视频hdxx| 99精国产麻豆久久婷婷| 九色亚洲精品在线播放| 美女脱内裤让男人舔精品视频| 免费高清在线观看日韩| 久久中文字幕一级| 国产欧美日韩一区二区精品| bbb黄色大片| 亚洲第一欧美日韩一区二区三区 | 中文字幕另类日韩欧美亚洲嫩草| 桃红色精品国产亚洲av| 欧美一级毛片孕妇| 一区二区三区四区激情视频| 午夜两性在线视频| 欧美精品高潮呻吟av久久| 亚洲精品美女久久av网站| 精品免费久久久久久久清纯 | 久久久欧美国产精品| 精品国产一区二区三区久久久樱花| 啦啦啦啦在线视频资源| 国产成人系列免费观看| 19禁男女啪啪无遮挡网站| 久久人妻熟女aⅴ| 视频在线观看一区二区三区| 在线观看免费午夜福利视频| 日本撒尿小便嘘嘘汇集6| 人人妻人人添人人爽欧美一区卜| 一区二区三区乱码不卡18| 老司机靠b影院| 亚洲精品美女久久av网站| 性高湖久久久久久久久免费观看| 免费人妻精品一区二区三区视频| 欧美精品av麻豆av| 岛国在线观看网站| 男女午夜视频在线观看| 在线观看一区二区三区激情| 国产一区二区激情短视频 | 精品久久久精品久久久| 曰老女人黄片| 久久狼人影院| 精品国产一区二区三区久久久樱花| 一区在线观看完整版| svipshipincom国产片| 亚洲av片天天在线观看| 国产精品免费视频内射| 国产精品麻豆人妻色哟哟久久| 欧美激情极品国产一区二区三区| a在线观看视频网站| 男女国产视频网站| 99国产精品一区二区三区| 性色av一级| 久久亚洲国产成人精品v| 天堂中文最新版在线下载| 如日韩欧美国产精品一区二区三区| 亚洲国产av影院在线观看| 成人影院久久| 又黄又粗又硬又大视频| 国产成人欧美| 日本精品一区二区三区蜜桃| 性色av乱码一区二区三区2| 欧美成狂野欧美在线观看| 日本猛色少妇xxxxx猛交久久| 欧美性长视频在线观看| 日本五十路高清| 男女床上黄色一级片免费看| 日韩有码中文字幕| 中文字幕最新亚洲高清| 精品久久久精品久久久| 一二三四在线观看免费中文在| 亚洲国产看品久久| 亚洲黑人精品在线| 中文精品一卡2卡3卡4更新| 在线观看舔阴道视频| 蜜桃国产av成人99| 欧美日韩中文字幕国产精品一区二区三区 | 欧美日韩亚洲综合一区二区三区_| 欧美日韩一级在线毛片| 精品高清国产在线一区| 日韩一区二区三区影片| 美女午夜性视频免费| 美女大奶头黄色视频| 日韩制服骚丝袜av| 色播在线永久视频| 妹子高潮喷水视频| 欧美变态另类bdsm刘玥| 亚洲七黄色美女视频| 午夜久久久在线观看| 99国产精品一区二区三区| 无遮挡黄片免费观看| 日本a在线网址| 女人精品久久久久毛片| 久久久久网色| 久久精品国产亚洲av香蕉五月 | 午夜免费成人在线视频| 黑丝袜美女国产一区| 精品熟女少妇八av免费久了| 国产成人免费无遮挡视频| 精品人妻1区二区| 亚洲av电影在线进入| 一区二区av电影网| 国产成人a∨麻豆精品| 国产精品免费视频内射| 日韩精品免费视频一区二区三区| 亚洲精品粉嫩美女一区| 亚洲人成电影免费在线| 国产成+人综合+亚洲专区| 日韩制服丝袜自拍偷拍| 国产一区二区激情短视频 | 亚洲 国产 在线| 男女之事视频高清在线观看| 亚洲欧洲日产国产| 欧美激情久久久久久爽电影 | 亚洲精品久久午夜乱码| 每晚都被弄得嗷嗷叫到高潮| 成在线人永久免费视频| 亚洲欧美成人综合另类久久久| 亚洲色图综合在线观看| 亚洲专区中文字幕在线| 人人妻人人澡人人爽人人夜夜| 亚洲av男天堂| 免费日韩欧美在线观看| 国产视频一区二区在线看| avwww免费| 一区二区三区精品91| 久久99热这里只频精品6学生| 十八禁高潮呻吟视频| 日韩视频一区二区在线观看| 一区在线观看完整版| 亚洲色图综合在线观看| 亚洲 国产 在线| 在线天堂中文资源库| 日本wwww免费看| 亚洲av成人一区二区三| 免费黄频网站在线观看国产| 精品免费久久久久久久清纯 | 一级毛片电影观看| 咕卡用的链子| 日本91视频免费播放| 国内毛片毛片毛片毛片毛片| 在线av久久热| 欧美人与性动交α欧美软件| 99久久人妻综合| 午夜福利在线免费观看网站| 99精品久久久久人妻精品| bbb黄色大片| 精品少妇内射三级| 热99国产精品久久久久久7| 久久精品熟女亚洲av麻豆精品| 久久免费观看电影| 日本av手机在线免费观看| 自线自在国产av| 久久人妻熟女aⅴ| 2018国产大陆天天弄谢| 久久精品亚洲熟妇少妇任你| 亚洲国产看品久久| 日本a在线网址| av网站在线播放免费| 欧美精品啪啪一区二区三区 | 国产精品av久久久久免费| 精品一区二区三区四区五区乱码| 精品高清国产在线一区| 日本91视频免费播放| 一级片免费观看大全| 午夜成年电影在线免费观看| 亚洲精品第二区| 亚洲成人国产一区在线观看| 在线观看免费午夜福利视频| 色播在线永久视频| 国产xxxxx性猛交| av超薄肉色丝袜交足视频| 亚洲情色 制服丝袜| 99热全是精品| 国产一区二区三区综合在线观看| 波多野结衣一区麻豆| 日日爽夜夜爽网站| 99久久人妻综合| 欧美激情高清一区二区三区| 国产亚洲欧美精品永久| 丝袜脚勾引网站| 亚洲,欧美精品.| 日韩制服骚丝袜av| 美女国产高潮福利片在线看| 午夜福利一区二区在线看| 制服诱惑二区| 咕卡用的链子| 欧美黑人精品巨大| 亚洲色图 男人天堂 中文字幕| 王馨瑶露胸无遮挡在线观看| 亚洲第一av免费看| 国产成人影院久久av| 精品一区二区三区av网在线观看 | 午夜91福利影院| 国产av一区二区精品久久| 人人妻人人澡人人看| 国产亚洲精品av在线| av国产免费在线观看| 午夜免费激情av| 麻豆一二三区av精品| av天堂在线播放| 女生性感内裤真人,穿戴方法视频| 99热6这里只有精品| 久久精品亚洲精品国产色婷小说| 久久天堂一区二区三区四区| 给我免费播放毛片高清在线观看| 99热这里只有是精品50| 成人av一区二区三区在线看| 国产一区二区激情短视频| 亚洲欧美日韩高清专用| 国产亚洲精品av在线| 少妇人妻一区二区三区视频| 久久天躁狠狠躁夜夜2o2o| 国产aⅴ精品一区二区三区波| 成人高潮视频无遮挡免费网站| 精品国产亚洲在线| 国产欧美日韩精品亚洲av| 亚洲七黄色美女视频| 国产精品 欧美亚洲| 欧美成人午夜精品| 18禁观看日本| 亚洲精品粉嫩美女一区| 亚洲成a人片在线一区二区| 曰老女人黄片| a级毛片a级免费在线| 成人特级黄色片久久久久久久| 久久中文字幕一级| 中文资源天堂在线| 欧美性猛交╳xxx乱大交人| 天堂av国产一区二区熟女人妻 | xxxwww97欧美| av天堂在线播放| 国内精品久久久久久久电影| 成人国产一区最新在线观看| 精品欧美国产一区二区三| 欧美绝顶高潮抽搐喷水| 黄频高清免费视频| 无限看片的www在线观看| 国产一区二区三区在线臀色熟女| 成人欧美大片| 亚洲精品美女久久久久99蜜臀| 50天的宝宝边吃奶边哭怎么回事| 在线免费观看的www视频| www.熟女人妻精品国产| 欧美黑人精品巨大| 校园春色视频在线观看| 少妇被粗大的猛进出69影院| 法律面前人人平等表现在哪些方面| 男女下面进入的视频免费午夜| 色综合欧美亚洲国产小说| 男人舔女人下体高潮全视频| 可以免费在线观看a视频的电影网站| 亚洲国产精品sss在线观看| 首页视频小说图片口味搜索| 午夜福利18| 国产v大片淫在线免费观看| 搡老岳熟女国产| 一卡2卡三卡四卡精品乱码亚洲| 亚洲自拍偷在线| 国产精品亚洲av一区麻豆| 成人三级黄色视频| 国产精品久久电影中文字幕| 久久婷婷人人爽人人干人人爱| 蜜桃久久精品国产亚洲av| 国产免费av片在线观看野外av| svipshipincom国产片| 亚洲成人中文字幕在线播放| 18禁黄网站禁片午夜丰满| 国内精品久久久久精免费| 很黄的视频免费| 国产精品一区二区精品视频观看| 床上黄色一级片| 99久久精品国产亚洲精品| 国产精品一区二区三区四区免费观看 | 国产97色在线日韩免费| 国产精品98久久久久久宅男小说| 麻豆国产97在线/欧美 | 两性午夜刺激爽爽歪歪视频在线观看 | 嫁个100分男人电影在线观看| 国产伦在线观看视频一区| 国产精品,欧美在线| 亚洲一区高清亚洲精品| 黄色 视频免费看| 亚洲五月天丁香| 国产亚洲精品综合一区在线观看 | 精品久久蜜臀av无| 妹子高潮喷水视频| 国产黄色小视频在线观看| 人妻丰满熟妇av一区二区三区| 亚洲美女黄片视频| a在线观看视频网站| 成人永久免费在线观看视频| 国产三级黄色录像| 欧美最黄视频在线播放免费| 女生性感内裤真人,穿戴方法视频| 中文字幕高清在线视频| 久久中文字幕一级| 国产精品国产高清国产av| 国产成人精品无人区| 精品不卡国产一区二区三区| 制服诱惑二区| 午夜免费成人在线视频| 久久草成人影院| 啦啦啦免费观看视频1| 午夜成年电影在线免费观看| 国产精品免费一区二区三区在线| 脱女人内裤的视频| 叶爱在线成人免费视频播放| 亚洲一区中文字幕在线| 久久精品成人免费网站| 成人av一区二区三区在线看| 亚洲精品在线美女| 欧美乱妇无乱码| 精品电影一区二区在线| 精品一区二区三区四区五区乱码| 久久久国产欧美日韩av| www.精华液| 一级a爱片免费观看的视频| 亚洲精品国产一区二区精华液| 99久久精品热视频| 国内久久婷婷六月综合欲色啪| 日韩欧美三级三区| 一区二区三区国产精品乱码| 成人午夜高清在线视频| 精品少妇一区二区三区视频日本电影| 欧美成狂野欧美在线观看| 波多野结衣高清无吗| 国产单亲对白刺激| 久久这里只有精品中国| 级片在线观看| 哪里可以看免费的av片| 亚洲国产日韩欧美精品在线观看 | 国产av不卡久久| 9191精品国产免费久久| 久久精品国产亚洲av高清一级| 一二三四在线观看免费中文在| 成年人黄色毛片网站| 成人国产一区最新在线观看| aaaaa片日本免费| 日韩成人在线观看一区二区三区| 在线播放国产精品三级| 床上黄色一级片| 搡老熟女国产l中国老女人| 欧美zozozo另类| 在线观看66精品国产| 成人手机av| 日韩欧美 国产精品| 一级作爱视频免费观看| 亚洲精品在线观看二区| 特大巨黑吊av在线直播| 免费在线观看完整版高清| 精品国产乱子伦一区二区三区| 国产av一区二区精品久久| 欧洲精品卡2卡3卡4卡5卡区| 88av欧美| 国产探花在线观看一区二区| 91成年电影在线观看| 亚洲自偷自拍图片 自拍| 香蕉丝袜av| 亚洲专区国产一区二区| 日韩有码中文字幕| 国产精品香港三级国产av潘金莲| 国产成人啪精品午夜网站| 变态另类丝袜制服| 欧美 亚洲 国产 日韩一| 久久精品成人免费网站| 久久久久久人人人人人| 午夜老司机福利片| av国产免费在线观看| 精品久久久久久久久久久久久| 免费高清视频大片| 日韩欧美在线乱码| 母亲3免费完整高清在线观看| 一级黄色大片毛片| 真人一进一出gif抽搐免费| 在线观看www视频免费| 香蕉国产在线看| 97碰自拍视频| 欧美中文日本在线观看视频| 丰满的人妻完整版| 波多野结衣巨乳人妻| 露出奶头的视频| 成人精品一区二区免费| 亚洲精品国产精品久久久不卡| 老司机午夜福利在线观看视频| 国产亚洲欧美98| 国产精品综合久久久久久久免费| 欧美在线一区亚洲| 欧美绝顶高潮抽搐喷水| 亚洲,欧美精品.| 欧美zozozo另类| 麻豆成人午夜福利视频| 久久久久性生活片| 无遮挡黄片免费观看| 精品午夜福利视频在线观看一区| 成人三级黄色视频| 夜夜夜夜夜久久久久| 伊人久久大香线蕉亚洲五| 19禁男女啪啪无遮挡网站| 国产久久久一区二区三区| 中文字幕人成人乱码亚洲影| 欧美另类亚洲清纯唯美| 国产精品99久久99久久久不卡| 亚洲五月天丁香| 老司机在亚洲福利影院| 欧美+亚洲+日韩+国产| 国产成人精品久久二区二区91| 亚洲中文日韩欧美视频| 国产精品久久久人人做人人爽| 免费看日本二区| 老司机靠b影院| 国产真实乱freesex| 999久久久精品免费观看国产| 日本 欧美在线| 国产精品1区2区在线观看.| 亚洲五月天丁香| 欧美激情久久久久久爽电影| 亚洲人成77777在线视频| 精品一区二区三区av网在线观看| 欧美午夜高清在线| 中文在线观看免费www的网站 | 久久久久久久午夜电影| 一区二区三区国产精品乱码| 麻豆国产97在线/欧美 | 欧美激情久久久久久爽电影| 天堂√8在线中文| 亚洲欧洲精品一区二区精品久久久| 亚洲自拍偷在线| 亚洲第一欧美日韩一区二区三区| 97人妻精品一区二区三区麻豆| 黄色片一级片一级黄色片| 欧美+亚洲+日韩+国产| 三级毛片av免费| 亚洲熟女毛片儿| 香蕉国产在线看| 999久久久国产精品视频| 国产视频内射| 大型黄色视频在线免费观看| 麻豆国产av国片精品| 身体一侧抽搐| 亚洲色图 男人天堂 中文字幕| 在线观看舔阴道视频| 一二三四社区在线视频社区8| 观看免费一级毛片| 亚洲国产精品sss在线观看| 久久香蕉激情| 人妻夜夜爽99麻豆av| 日本一区二区免费在线视频| 美女免费视频网站| 国产97色在线日韩免费| 国产精品久久久av美女十八| 国产精品久久电影中文字幕| 男女午夜视频在线观看| 午夜激情福利司机影院| 18禁观看日本| 在线看三级毛片| 国产精品 欧美亚洲| 久久99热这里只有精品18| 国内久久婷婷六月综合欲色啪| 久久人人精品亚洲av| 国产视频内射| 五月伊人婷婷丁香| 久久九九热精品免费| 精品不卡国产一区二区三区| 欧美成狂野欧美在线观看| 男人舔奶头视频| 欧美色视频一区免费| 国产精品亚洲一级av第二区| 一边摸一边抽搐一进一小说| 麻豆成人av在线观看| 久久天躁狠狠躁夜夜2o2o| av福利片在线观看| 一本一本综合久久| 日韩免费av在线播放| 亚洲人成电影免费在线| 久久婷婷成人综合色麻豆| 成年人黄色毛片网站| 神马国产精品三级电影在线观看 | 亚洲avbb在线观看| 精品人妻1区二区|