• <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
    色在线成人网| 交换朋友夫妻互换小说| 日本vs欧美在线观看视频| 久久人人爽av亚洲精品天堂| 777久久人妻少妇嫩草av网站| 老司机靠b影院| 一区二区日韩欧美中文字幕| 在线看a的网站| 国产麻豆69| 中文字幕av电影在线播放| 日韩中文字幕欧美一区二区| 午夜福利,免费看| 亚洲五月天丁香| 国产aⅴ精品一区二区三区波| 丝袜美足系列| 亚洲全国av大片| 一二三四社区在线视频社区8| 亚洲人成电影观看| 丰满迷人的少妇在线观看| 亚洲国产精品合色在线| 欧美黑人精品巨大| 最近最新中文字幕大全免费视频| 熟女少妇亚洲综合色aaa.| 老司机深夜福利视频在线观看| 日本五十路高清| 日本精品一区二区三区蜜桃| 麻豆国产av国片精品| 国产激情久久老熟女| 香蕉久久夜色| 村上凉子中文字幕在线| cao死你这个sao货| 黄色丝袜av网址大全| 曰老女人黄片| 一a级毛片在线观看| 精品国产一区二区三区久久久樱花| 99久久国产精品久久久| 精品久久久精品久久久| 国产野战对白在线观看| а√天堂www在线а√下载 | 侵犯人妻中文字幕一二三四区| 日韩 欧美 亚洲 中文字幕| 国产av精品麻豆| 三级毛片av免费| 热re99久久国产66热| 91字幕亚洲| 亚洲精品国产精品久久久不卡| 视频在线观看一区二区三区| xxx96com| 大型黄色视频在线免费观看| 美女高潮喷水抽搐中文字幕| 99精品久久久久人妻精品| 亚洲在线自拍视频| 亚洲专区中文字幕在线| 久久中文字幕人妻熟女| 精品午夜福利视频在线观看一区| 国产乱人伦免费视频| 国产不卡一卡二| 最近最新中文字幕大全免费视频| 黑人猛操日本美女一级片| 成年人午夜在线观看视频| 免费少妇av软件| 51午夜福利影视在线观看| 久久精品国产99精品国产亚洲性色 | 美女视频免费永久观看网站| 日韩熟女老妇一区二区性免费视频| 别揉我奶头~嗯~啊~动态视频| 脱女人内裤的视频| 999精品在线视频| 99在线人妻在线中文字幕 | 欧美精品高潮呻吟av久久| 免费观看人在逋| 在线av久久热| 777久久人妻少妇嫩草av网站| 欧美黑人精品巨大| 中亚洲国语对白在线视频| 国产精品 欧美亚洲| 亚洲精品一卡2卡三卡4卡5卡| 久热这里只有精品99| 18禁美女被吸乳视频| 久久久久久亚洲精品国产蜜桃av| 欧美日韩精品网址| 亚洲男人天堂网一区| 欧美乱妇无乱码| 免费在线观看日本一区| 欧美+亚洲+日韩+国产| 大香蕉久久网| 女人被狂操c到高潮| 免费在线观看视频国产中文字幕亚洲| 久久久精品免费免费高清| 久久久久国内视频| 久久青草综合色| 日韩免费av在线播放| 国产精品一区二区在线观看99| e午夜精品久久久久久久| 国产成+人综合+亚洲专区| 久久久国产欧美日韩av| 1024视频免费在线观看| 丁香六月欧美| 国产精品久久久av美女十八| 又大又爽又粗| 黄色a级毛片大全视频| 丝袜在线中文字幕| 国产精品一区二区在线观看99| 久久久久久亚洲精品国产蜜桃av| 亚洲精品成人av观看孕妇| 亚洲男人天堂网一区| 欧美色视频一区免费| 国产精品电影一区二区三区 | 人妻久久中文字幕网| 国产精品99久久99久久久不卡| 午夜福利一区二区在线看| av网站免费在线观看视频| 黄色丝袜av网址大全| 动漫黄色视频在线观看| 大陆偷拍与自拍| 老司机午夜十八禁免费视频| 国产精品乱码一区二三区的特点 | 色在线成人网| 久久久久久久久久久久大奶| 亚洲在线自拍视频| 村上凉子中文字幕在线| 久久精品亚洲av国产电影网| 免费一级毛片在线播放高清视频 | 免费日韩欧美在线观看| 美女高潮喷水抽搐中文字幕| 国产精品亚洲av一区麻豆| 国产免费男女视频| 80岁老熟妇乱子伦牲交| 操出白浆在线播放| 欧美精品啪啪一区二区三区| 亚洲久久久国产精品| 麻豆国产av国片精品| 色播在线永久视频| 国产免费av片在线观看野外av| 久久九九热精品免费| 一级毛片高清免费大全| 国产又色又爽无遮挡免费看| 亚洲熟女精品中文字幕| 老司机深夜福利视频在线观看| 久久久精品免费免费高清| 如日韩欧美国产精品一区二区三区| 国产精品电影一区二区三区 | 国产片内射在线| 国产99白浆流出| 狂野欧美激情性xxxx| 一级a爱视频在线免费观看| 宅男免费午夜| 亚洲免费av在线视频| 下体分泌物呈黄色| 人人妻人人爽人人添夜夜欢视频| 精品无人区乱码1区二区| 成年人午夜在线观看视频| 少妇被粗大的猛进出69影院| 亚洲欧美激情综合另类| 在线观看免费视频日本深夜| 大陆偷拍与自拍| √禁漫天堂资源中文www| 一边摸一边做爽爽视频免费| 叶爱在线成人免费视频播放| 精品久久久久久久久久免费视频 | 国产精品香港三级国产av潘金莲| 欧美精品一区二区免费开放| 热99国产精品久久久久久7| 一进一出好大好爽视频| 女人精品久久久久毛片| 在线观看免费视频日本深夜| 91av网站免费观看| 久久青草综合色| 欧美av亚洲av综合av国产av| 波多野结衣一区麻豆| 91在线观看av| 精品乱码久久久久久99久播| 亚洲第一欧美日韩一区二区三区| 免费在线观看影片大全网站| 电影成人av| 看免费av毛片| 国产不卡一卡二| 欧美色视频一区免费| 日韩免费高清中文字幕av| 亚洲片人在线观看| 国产麻豆69| 国产精品一区二区精品视频观看| 亚洲熟妇熟女久久| 国产激情久久老熟女| 飞空精品影院首页| 老鸭窝网址在线观看| xxxhd国产人妻xxx| 黄色女人牲交| 国产精品亚洲av一区麻豆| 亚洲熟女毛片儿| 成人亚洲精品一区在线观看| 亚洲情色 制服丝袜| 日韩欧美国产一区二区入口| 十八禁人妻一区二区| 一级毛片高清免费大全| 免费在线观看完整版高清| 国产精品秋霞免费鲁丝片| 欧美不卡视频在线免费观看 | 女警被强在线播放| 9热在线视频观看99| 久久香蕉国产精品| 校园春色视频在线观看| av天堂在线播放| 久久精品成人免费网站| 久久精品国产综合久久久| 国产不卡av网站在线观看| 国产精品秋霞免费鲁丝片| 亚洲精品一卡2卡三卡4卡5卡| 五月开心婷婷网| www.精华液| ponron亚洲| 亚洲av欧美aⅴ国产| 日日摸夜夜添夜夜添小说| 美女午夜性视频免费| 两性夫妻黄色片| 99国产精品一区二区三区| 亚洲成av片中文字幕在线观看| 亚洲免费av在线视频| 久久久久久免费高清国产稀缺| 久久人妻福利社区极品人妻图片| 亚洲 欧美一区二区三区| 巨乳人妻的诱惑在线观看| 波多野结衣一区麻豆| 亚洲国产精品合色在线| 99精国产麻豆久久婷婷| 中文字幕av电影在线播放| 国产主播在线观看一区二区| 国产成人欧美在线观看 | 国产区一区二久久| ponron亚洲| 亚洲aⅴ乱码一区二区在线播放 | 日韩成人在线观看一区二区三区| 18禁黄网站禁片午夜丰满| 97人妻天天添夜夜摸| 女性生殖器流出的白浆| 人妻丰满熟妇av一区二区三区 | 免费在线观看亚洲国产| 999久久久精品免费观看国产| 日韩欧美免费精品| 91精品国产国语对白视频| 久热这里只有精品99| aaaaa片日本免费| 激情视频va一区二区三区| 欧美激情高清一区二区三区| 国产1区2区3区精品| 精品人妻在线不人妻| 免费在线观看亚洲国产| 午夜精品国产一区二区电影| 中文欧美无线码| www日本在线高清视频| 久久久水蜜桃国产精品网| 国产亚洲欧美精品永久| 在线视频色国产色| 亚洲欧美激情综合另类| 大型黄色视频在线免费观看| 久9热在线精品视频| 欧美激情久久久久久爽电影 | 免费在线观看完整版高清| 亚洲第一欧美日韩一区二区三区| 国产日韩欧美亚洲二区| 欧美乱妇无乱码| 手机成人av网站| 亚洲性夜色夜夜综合| 亚洲在线自拍视频| 黄片小视频在线播放| 一级片'在线观看视频| 成人国产一区最新在线观看| 国产av又大| 老司机在亚洲福利影院| 午夜亚洲福利在线播放| 欧美日韩国产mv在线观看视频| 99热只有精品国产| www.熟女人妻精品国产| 午夜福利欧美成人| 女人精品久久久久毛片| 99国产极品粉嫩在线观看| 久久九九热精品免费| 欧美日韩精品网址| 69av精品久久久久久| 看片在线看免费视频| 国产一卡二卡三卡精品| 免费av中文字幕在线| 亚洲欧美色中文字幕在线| 欧美日韩一级在线毛片| 999久久久精品免费观看国产| 欧美精品人与动牲交sv欧美| 亚洲一区二区三区不卡视频| 黑丝袜美女国产一区| 亚洲欧美日韩高清在线视频| 免费在线观看完整版高清| 国产男女内射视频| 啦啦啦免费观看视频1| 性少妇av在线| 亚洲精品在线观看二区| 极品人妻少妇av视频| 黄色 视频免费看| 日韩大码丰满熟妇| 女同久久另类99精品国产91| 一级,二级,三级黄色视频| 国产免费男女视频| 午夜亚洲福利在线播放| 一本大道久久a久久精品| 一二三四在线观看免费中文在| 亚洲av电影在线进入| 夜夜躁狠狠躁天天躁| 黄色女人牲交| 高清在线国产一区| 一a级毛片在线观看| 十分钟在线观看高清视频www| 三上悠亚av全集在线观看| 国产精品成人在线| 一区二区三区激情视频| 国产成人影院久久av| 国产精品乱码一区二三区的特点 | 久久香蕉国产精品| 中文字幕av电影在线播放| 一本一本久久a久久精品综合妖精| 精品欧美一区二区三区在线| 波多野结衣av一区二区av| 免费在线观看亚洲国产| 午夜免费成人在线视频| 夜夜夜夜夜久久久久| 嫩草影视91久久| 男人操女人黄网站| 欧美 日韩 精品 国产| 欧美成狂野欧美在线观看| 女同久久另类99精品国产91| 国产精品久久视频播放| 亚洲精品一卡2卡三卡4卡5卡| 亚洲一码二码三码区别大吗| 免费日韩欧美在线观看| 免费女性裸体啪啪无遮挡网站| 看黄色毛片网站| 亚洲精品美女久久av网站| 国产在线精品亚洲第一网站| 在线永久观看黄色视频| 国产高清videossex| 精品国产一区二区三区四区第35| 欧美日韩亚洲国产一区二区在线观看 | 久久精品人人爽人人爽视色| 午夜日韩欧美国产| 国产三级黄色录像| 看免费av毛片| 日韩熟女老妇一区二区性免费视频| 亚洲色图综合在线观看| 国产午夜精品久久久久久| 国产成人av激情在线播放| 精品乱码久久久久久99久播| 我的亚洲天堂| 亚洲精品国产精品久久久不卡| 精品少妇久久久久久888优播| 国产91精品成人一区二区三区| 在线播放国产精品三级| 在线免费观看的www视频| 成人18禁在线播放| 欧美激情 高清一区二区三区| 色尼玛亚洲综合影院| av在线播放免费不卡| 欧美一级毛片孕妇| 欧美黑人精品巨大| 制服人妻中文乱码| 99re6热这里在线精品视频| 国产精品一区二区免费欧美| 99re6热这里在线精品视频| 午夜免费鲁丝| 成年人免费黄色播放视频| 精品国产乱子伦一区二区三区| 一本一本久久a久久精品综合妖精| 最近最新中文字幕大全电影3 | 一级片'在线观看视频| 国产一区在线观看成人免费| 欧美国产精品一级二级三级| 精品国产美女av久久久久小说| 成熟少妇高潮喷水视频| 久久精品aⅴ一区二区三区四区| 婷婷丁香在线五月| av欧美777| 法律面前人人平等表现在哪些方面| 日韩 欧美 亚洲 中文字幕| av网站在线播放免费| 99re在线观看精品视频| 高清在线国产一区| 天堂动漫精品| 日日摸夜夜添夜夜添小说| 国产精品.久久久| 久久 成人 亚洲| 人人妻人人澡人人爽人人夜夜| 一级毛片精品| 99国产综合亚洲精品| 日韩免费高清中文字幕av| 性色av乱码一区二区三区2| 欧美 亚洲 国产 日韩一| 精品免费久久久久久久清纯 | tocl精华| 男人舔女人的私密视频| 老司机午夜十八禁免费视频| 丁香欧美五月| 法律面前人人平等表现在哪些方面| 80岁老熟妇乱子伦牲交| 黄片播放在线免费| 成人免费观看视频高清| 欧美乱色亚洲激情| 母亲3免费完整高清在线观看| 久热爱精品视频在线9| 国产黄色免费在线视频| 老司机靠b影院| 黄色片一级片一级黄色片| 嫩草影视91久久| 国产成人免费观看mmmm| 欧美在线黄色| 色尼玛亚洲综合影院| 黄片小视频在线播放| 一级毛片女人18水好多| 日韩视频一区二区在线观看| 久久精品亚洲熟妇少妇任你| 最近最新免费中文字幕在线| 国产一区二区三区视频了| 亚洲一卡2卡3卡4卡5卡精品中文| 乱人伦中国视频| 久久久久视频综合| 他把我摸到了高潮在线观看| 91老司机精品| 午夜福利在线免费观看网站| 国产精品久久久久成人av| 国产成人精品无人区| 亚洲色图av天堂| 午夜福利在线观看吧| 亚洲熟女毛片儿| 99香蕉大伊视频| 婷婷精品国产亚洲av在线 | 亚洲,欧美精品.| 欧美最黄视频在线播放免费 | 国产精品久久视频播放| 美女视频免费永久观看网站| aaaaa片日本免费| 国产xxxxx性猛交| √禁漫天堂资源中文www| 欧美老熟妇乱子伦牲交| 在线观看免费高清a一片| 久久中文字幕人妻熟女| 亚洲国产精品合色在线| 国产精品偷伦视频观看了| 一进一出抽搐gif免费好疼 | 高清av免费在线| 精品人妻在线不人妻| 亚洲精品乱久久久久久| 涩涩av久久男人的天堂| 操出白浆在线播放| 在线观看一区二区三区激情| av天堂久久9| 国产欧美日韩精品亚洲av| 99久久99久久久精品蜜桃| 亚洲国产精品一区二区三区在线| 亚洲成人免费电影在线观看| 午夜老司机福利片| av在线播放免费不卡| 欧美久久黑人一区二区| 夜夜爽天天搞| 9191精品国产免费久久| 国产精品久久久久久人妻精品电影| 国产亚洲av高清不卡| 欧美日韩一级在线毛片| 国产成人系列免费观看| 最近最新中文字幕大全电影3 | 最新美女视频免费是黄的| 欧美激情极品国产一区二区三区| 亚洲自偷自拍图片 自拍| 麻豆av在线久日| 少妇粗大呻吟视频| 久久久国产成人精品二区 | 色精品久久人妻99蜜桃| 久久婷婷成人综合色麻豆| 美女 人体艺术 gogo| 国产又爽黄色视频| 精品久久久久久久久久免费视频 | 国产精品 国内视频| 不卡一级毛片| 久久久国产成人精品二区 | 大型av网站在线播放| 黄网站色视频无遮挡免费观看| 一进一出抽搐动态| avwww免费| 黄色怎么调成土黄色| 99国产综合亚洲精品| 国产精品一区二区精品视频观看| av福利片在线| 丁香欧美五月| 美国免费a级毛片| 电影成人av| 成人av一区二区三区在线看| 在线av久久热| 在线免费观看的www视频| 在线国产一区二区在线| 人成视频在线观看免费观看| 999久久久精品免费观看国产| 免费不卡黄色视频| 人妻 亚洲 视频| 丁香六月欧美| 交换朋友夫妻互换小说| 午夜视频精品福利| 成熟少妇高潮喷水视频| 欧美日本中文国产一区发布| www日本在线高清视频| 成人手机av| 极品少妇高潮喷水抽搐| 国产精品永久免费网站| 美女高潮到喷水免费观看| 欧美精品啪啪一区二区三区| tocl精华| 免费少妇av软件| 精品一区二区三区av网在线观看| 欧美激情高清一区二区三区| 激情视频va一区二区三区| av片东京热男人的天堂| a级片在线免费高清观看视频| 国产成人精品久久二区二区免费| 日韩人妻精品一区2区三区| 每晚都被弄得嗷嗷叫到高潮| 老司机午夜福利在线观看视频| 99热只有精品国产| 又黄又粗又硬又大视频| 欧美精品一区二区免费开放| 亚洲国产毛片av蜜桃av| 午夜免费鲁丝| 男女免费视频国产| 婷婷成人精品国产| 人妻久久中文字幕网| 天堂√8在线中文| 欧美精品啪啪一区二区三区| 久久久国产成人免费| 精品一区二区三卡| 啦啦啦在线免费观看视频4| 亚洲av第一区精品v没综合| 婷婷成人精品国产| av网站在线播放免费| 19禁男女啪啪无遮挡网站| 国产精品国产高清国产av | 国产真人三级小视频在线观看| 午夜免费鲁丝| 男女下面插进去视频免费观看| 12—13女人毛片做爰片一| 大片电影免费在线观看免费| 国产三级黄色录像| 一二三四在线观看免费中文在| 午夜老司机福利片| 国产不卡av网站在线观看| 搡老熟女国产l中国老女人| 亚洲五月婷婷丁香| 午夜福利视频在线观看免费| 国产成人av激情在线播放| 日本a在线网址| 每晚都被弄得嗷嗷叫到高潮| 国产亚洲欧美精品永久| 女性生殖器流出的白浆| 亚洲欧美日韩另类电影网站| 色婷婷av一区二区三区视频| 一本一本久久a久久精品综合妖精| 亚洲自偷自拍图片 自拍| 麻豆国产av国片精品| 日韩精品免费视频一区二区三区| 啦啦啦在线免费观看视频4| 一级片免费观看大全| 国产深夜福利视频在线观看| 久久ye,这里只有精品| 国内毛片毛片毛片毛片毛片| 国产一区二区三区视频了| 人妻久久中文字幕网| 国产精品.久久久| 在线观看66精品国产| 黄色片一级片一级黄色片| 欧美午夜高清在线| 中文字幕精品免费在线观看视频| 久9热在线精品视频| 亚洲片人在线观看| 亚洲欧美一区二区三区久久| 欧美成狂野欧美在线观看| 国产不卡一卡二| 国产成人精品久久二区二区免费| 高清毛片免费观看视频网站 | 高清视频免费观看一区二区| 女性被躁到高潮视频| 三级毛片av免费| 国产精品免费一区二区三区在线 | 午夜福利免费观看在线| 国产亚洲精品久久久久久毛片 | 亚洲男人天堂网一区| 又黄又爽又免费观看的视频| 久久久精品国产亚洲av高清涩受| 亚洲第一青青草原| 岛国在线观看网站| 精品一区二区三区四区五区乱码| 亚洲,欧美精品.| 岛国在线观看网站| 免费在线观看亚洲国产| 亚洲成av片中文字幕在线观看| 9热在线视频观看99| 色婷婷av一区二区三区视频| www.999成人在线观看| 水蜜桃什么品种好| 久久久久久亚洲精品国产蜜桃av| 亚洲精品美女久久av网站| 在线观看66精品国产| 久久精品国产99精品国产亚洲性色 | 中文字幕精品免费在线观看视频| 欧美日韩一级在线毛片| 午夜免费观看网址| 岛国在线观看网站| 亚洲欧美激情综合另类| 亚洲精品美女久久久久99蜜臀| 精品福利观看| 操美女的视频在线观看| 亚洲精品美女久久av网站| 精品少妇一区二区三区视频日本电影| 老司机午夜福利在线观看视频| 视频区欧美日本亚洲|