999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

一種改進的虛擬手術中人體軟組織形變方法

2015-10-12 05:22:56何巍王魏平師為禮苗語何飛楊華民
關鍵詞:結構手術模型

何巍,王魏平,師為禮,苗語,何飛,楊華民

(長春理工大學 計算機科學與技術學院,長春 130022)

一種改進的虛擬手術中人體軟組織形變方法

何巍,王魏平,師為禮,苗語,何飛,楊華民

(長春理工大學計算機科學與技術學院,長春130022)

虛擬手術是將計算機虛擬現(xiàn)實技術應用于現(xiàn)代醫(yī)學,利用醫(yī)學圖像數(shù)據(jù),重構出虛擬人體軟組織模型,利用觸覺交互設備,在視覺和觸覺上為使用者提供手術場景的模擬和仿真。其中,軟組織形變的仿真和建模是虛擬手術系統(tǒng)核心的研究內容之一。針對虛擬手術中軟組織形變仿真,提出了一種改進的彈簧—質點模型。在正四邊形網(wǎng)格拓撲結構表面模型中,增加了結構彈簧、剪切彈簧和彎曲彈簧,使表面模型具有體模型的特征。該模型繼承了傳統(tǒng)彈簧—質點模型的建??臁⒃砗唵?、仿真速度快等優(yōu)點,同時還具有控制形變區(qū)域的能力,從而提高了仿真模型的精確度,滿足交互系統(tǒng)的真實性和穩(wěn)定性需求。實驗表明該方法有效增強了形變仿真的體積特征,適合于軟組織形變較大時的實時仿真,滿足了虛擬手術系統(tǒng)實時性的要求。

虛擬現(xiàn)實;彈簧—質點模型;軟組織形變;碰撞檢測

隨著計算機科學技術的不斷發(fā)展,虛擬現(xiàn)實技術(Virtual Reality Technology VRT)已被廣泛應用于現(xiàn)代科學研究的各個領域[1,2],其中利用虛擬現(xiàn)實技術進行生物醫(yī)學仿真是當前的一個研究熱點,虛擬手術系統(tǒng)是其中的典型應用。

虛擬手術系統(tǒng)是對醫(yī)學數(shù)據(jù)進行可視化并實現(xiàn)實時交互,模擬手術過程中組織器官變形、感官反饋等可能遇到的各種現(xiàn)象的虛擬現(xiàn)實應用系統(tǒng)。其中在軟組織器官形變方面,生物力學領域的研究提出了很多復雜的數(shù)學模型,能夠比較準確地描述軟組織的形變,軟組織大多數(shù)物理模型都是建立在彈性力學理論上的。Santanu等人的基于有限元方法(Finite Element Model,F(xiàn)EM)提出了具有生物特征的物理形變模型,雖然精確度很高,但是網(wǎng)格數(shù)量很多,導致模型實時性差[3];Zhong等人提出具有物理化學特征的表面模型,但很難進行復雜的虛擬手術操作[4];陳衛(wèi)東等改進了基于表面正六邊形模型,雖然計算量小,實時性高,但是仿真精度有待提高[5];王瑞藝提出的基于無網(wǎng)格SPH的物理模型,雖然能夠動態(tài)地反映軟組織的實時形變,但交互方面需要進一步改善[6];李艷東,朱玲等人結合ALMSDM的特點提出了頂點法向量局部更新與預計算策略,提高了實時性,但在大范圍變形的情況下,穩(wěn)定性不足[7];Wang等用邊界元法建立模型,雖對模型的邊界進行了離散,簡化了計算,但穩(wěn)定性較差[8]。因此,在保證模型形變精確度的同時,提高交互的實時性,是我們虛擬手術操作中需要解決的首要問題。傳統(tǒng)的基于表面網(wǎng)格拓撲結構的物理模型,存在兩方面問題,沒有體模型的特征或者形變真實性不足。

本論文在保證實時性的同時,在虛擬的人體軟組織表面模型中設置了結構彈簧、剪切彈簧和彎曲彈簧,使模型具備了體模型的特征,從而提高了仿真精度,達到了更好的仿真效果。同時保證形變實時性高,很好地滿足虛擬手術中交互過程的真實性、實時性和穩(wěn)定性需求。

1 軟組織形變物理模型

彈簧-質點模型(Mass Spring Method,MSM)是一種經(jīng)典的物理建模方法,其主要思想是把仿真的軟組織用質點離散,質點與質點之間用滿足胡克定律的彈簧連接。該模型通常基于一定的幾何拓撲結構,形變過程是當一個質點受外力作用時,產(chǎn)生的應力作用在其它相鄰的質點上,同時通過彈簧傳遞力,帶動其相鄰的質點運動。模型的建立比較簡單,形變的計算量較低,能夠很好滿足實時性的要求。

1.1軟組織的生物力學特征

對人體軟組織器官(如血管、肝臟、皮膚、氣管等)形變特性的研究屬于生物力學范疇[9]。隨著生物力學研究的發(fā)展,軟組織的生物力學特征通常表現(xiàn)為非線性、粘彈性、不均勻性和各向異性等[10]。因此,虛擬手術仿真系統(tǒng)中的軟組織建模必須考慮上述生物力學特性,并根據(jù)實際情況做相應的簡化處理。通常采用彈性模量、阻尼系數(shù)等物理量來表示軟組織的各種生物力學特征。

1.2彈簧—質點模型建立

假設彈簧—質點模型由n個質點組成,在任意時刻t每個質點Ni的運動方程都滿足牛頓第二定律,因此每個質點Ni都具有如下的動力學微分方程(1):

其中,mi表示第i個質點質量;xi表示質點Ni的形變位移;ci表示阻尼系數(shù);kij表示彈簧的彈性系數(shù);lij表示質點Ni和Nj間彈簧的初始長度,lij0表示形變后的長度;Fext表示質點Ni所受的外力的綜合。

在彈簧—質點模型模擬人體軟組織形變時,力的作用導致質點發(fā)生位移。力的構成,除了外力之外,還有彈簧的變形產(chǎn)生的內力。不同的彈簧布局在仿真過程中會產(chǎn)生不同的力,質點也會產(chǎn)生不同的位移,使模型產(chǎn)生不同的變形效果。當采用彈簧—質點模型進行軟組織變形仿真實驗時,核心之一是建立合適的幾何拓撲結構。

本文對人體的皮膚進行仿真實驗,采用面模型。因為人體的皮膚質量分布相對均勻,其網(wǎng)格可均勻劃分。網(wǎng)格劃分的密度應根據(jù)仿真需求的精度和交互的實時性來確定。網(wǎng)格的拓撲結構從不同的精度考慮采用圖1(b)、(c)、(d)三種方式[11]。

圖1 人體軟組織面模型拓撲結構

為提高仿真精度,我們選用圖1(d)中拓撲結構作為物理模型。彈簧可分為3類:結構彈簧、剪切彈簧和彎曲彈簧,如圖2所示。

圖2 彈簧—質點模型的三種彈簧類型

(1)結構彈簧:質點Ni,j和Ni+1,j之間和質點Ni,j和Ni,j+1之間的彈簧,是防止質點之間過度的拉壓形變。

(2)剪切彈簧:質點Ni,j和Ni+1,j+1之間和質點Ni,j+1和 Ni+1,j之間的彈簧,是阻止結構的斜向過度變形,給結構一個剪切剛度。

(3)彎曲彈簧:質點Ni,j和Ni+2,j之間和質點Ni,j和Ni,j+2之間的彈簧,是為了防止結構過度的不自然彎曲。

2 軟組織模型算法分析

2.1動力學方程

基于上述彈簧—質點模型拓撲結構,結合受外力情況下模型的運動規(guī)律,對質點建立動力學方程,通過求方程的解來實現(xiàn)軟組織的形變。

力學模型的建立基于以下兩點:(1)模型能夠達到實時的仿真效果;(2)在現(xiàn)有的硬件條件下能夠實現(xiàn)相關算法的復雜度要求。現(xiàn)有的基于物理模型的仿真算法一般采用兩種基本方法:力法和能量法。為了滿足仿真的實時交互需求,本文采用力法。以下是對彈簧—質點模型的動力學分析[12,13]:

由牛頓第二定律可得,模型中每個質點Ni的運動微分方程由它所受到的合力(內力和外力組成)決定,對于質點N 有:

其中,F(xiàn)ext(i)是質點Ni所受的外力,F(xiàn)int(i)是質點Ni所受的內力(彈性力、阻尼力等)。根據(jù)胡克定律,質點Ni所受彈簧的拉(壓)力如公式(3)所示:

其中,Kijs是質點Ni和Nj間的彈性系數(shù),Iij是該時刻彈簧的矢量,Xi和Xj是質點Ni和Nj的位置。

質點Ni和Nj之間的彈簧阻尼力如下:

其中,Kijd是質點Ni和Nj之間的彈簧阻尼系數(shù)。

根據(jù)公式(3)和(4),可以得到質點Ni所受的內力:

其中,K是與質點Ni連接的質點總數(shù)。

一個由N個質點組成的模型,所受的內力由結構力,彎曲力和剪切力組成;可以表述為:

其中,F(xiàn)struct、Fshearing、Fbending分別是與質點Ni連接的彈簧產(chǎn)生的結構力、剪切力和彎曲力,k1、k2、k3分別表示與質點Ni相連的結構彈簧、剪切彈簧和彎曲彈簧的數(shù)量。

一般情況下,質點Ni所受到的外力包括重力和拉力,則受到的外力可以表述為:

對于上述的動力學微分方程,模型x并不是現(xiàn)實給出的,為了方便求解,需要對二階微分方程進行處理,將其轉化為一階常微分方程進行求解。則,每個質點Ni的速度和位移都可以用以下公式表示:

2.2模型的數(shù)值分析

首先考慮一階常微分方程的初值求解問題:

其中f( )x,y為已知函數(shù);y0為初值。

數(shù)值積分的原理就是,在區(qū)間[a,b]上取n+1個節(jié) a=x0<x1<x2<…<xn=b。 hi=xi-xi-1(1,2,…) 表示相鄰兩個節(jié)點的間距,稱為步長。這些hi可以不相等,我們取相等考慮。則將這些節(jié)點離散化,就可以將初值問題轉化成離散變量的相關問題。把相應問題的解yn作為y(xn)的近似值。這樣解得yn就是上述初值問題在節(jié)點xn的數(shù)值解,不同的離散方法產(chǎn)生不同的結果。

一階常微分方程初值問題的幾種常用的數(shù)值積分算法包括歐拉方法、中點方法以及四階龍格—庫塔方法等。對于質點數(shù)一定的軟組織模型,形變的實時性和穩(wěn)定性取決于算法。主要從以下幾個方面來考慮:

(1)步長:表示為了達到一定數(shù)值穩(wěn)定性所需要的時間離散化程度,表明方法的數(shù)值穩(wěn)定性;

(2)迭代次數(shù):反映了采用積分方法的計算量,影響總的計算速度;

(3)誤差:數(shù)值求解精度受誤差影響,包括截斷誤差、舍入誤差等;當確定步長后,微分方程的階次越高,截斷誤差越?。槐?列出了以上三種數(shù)值算法的截斷誤差。

表1 三種數(shù)值算法的截斷誤差

從計算精度上來看,步長h確定的情況下,四階龍格—庫塔法的計算精度最高,但是執(zhí)行一步所需的計算量也最大,影響了仿真形變的實時性。只從單步積分來看,步長設置越小,模型離散化程度就越高,局部截斷誤差也就越小。但是離散化程度越大,同樣意味著在一定范圍內需要的步數(shù)越多,不僅增加了計算量,還有可能引起舍入誤差的積累過大,導致精度變低。

3 碰撞檢測

碰撞檢測稱為接觸檢測或者干涉檢測,基本目的是確定兩個或兩個以上物體彼此之間是否發(fā)生接觸或穿透并計算出碰撞發(fā)生的位置。在虛擬手術中,需要碰撞檢測的場景主要有手術器械與人體器官之間、人體器官的自碰撞和手術器械的自碰撞。這些是整個虛擬手術系統(tǒng)中最重要的部分之一,決定著下一步的軟組織變形、切割、縫合等的實現(xiàn)。

如圖3所示,質點t0時刻在軟組織外部,而在下一個時間步長t0+Δt時刻已經(jīng)在軟組織內部,則質點和軟組織在t0+Δt和t0+Δt之間的時刻發(fā)生了碰撞。

圖3 碰撞對象

目前常見的碰撞檢測算法主要有:層次包圍盒法[14]和空間分解法[15]。

層次包圍盒的核心思想是用體積略大而幾何特征簡單的包圍盒近似表述復雜的幾何對象,只對包圍盒相交時包裹的對象進行相交實驗。此外,通過引入樹狀層次結構快速刪除不發(fā)生碰撞的部分,減少了不必要的相交實驗,提高檢測效率。層次包圍盒法應用較為廣泛,適用于復雜環(huán)境下的碰撞檢測。根據(jù)所采用的包圍體類型的不同可對層次包圍盒進行區(qū)分,比較典型的有沿坐標軸的包圍盒AABB、方向包圍盒OBB和包圍球等。

空間分解法的核心思想是將全部空間分解成體積相等的小單元格,只對同一單元格或相鄰單元格的幾何對象進行相交測試。比較典型的有八叉樹、BSP樹等。

虛擬手術仿真系統(tǒng)需要很好地滿足交互實時性的要求,因此,我們使用AABB層次包圍盒的改進算法作為本文虛擬手術仿真系統(tǒng)的碰撞檢測算法。

AABB方法[16]沿坐標軸的包圍盒是一種使用廣泛的碰撞檢測算法,包含幾何對象且邊平行于坐標軸的最小六面體。因此只需要六個標量就能描述一個AABB方法。

AABB樹是基于AABB方法自底至上動態(tài)更新生成的層次結構二叉樹。AABB樹算法的核心是通過遍歷兩個碰撞物體的AABB樹來判斷該位置是否發(fā)生碰撞。步驟如下:

步驟1:構造碰撞剛體和軟組織的AABB樹;

步驟2:檢測最大包圍盒是否相交;

步驟3:遞歸處理剛體和軟組織的AABB樹,如果發(fā)現(xiàn)子樹沒有發(fā)生相交,停止并得出結論沒有發(fā)生碰撞。如果發(fā)現(xiàn)子樹相交,則進一步處理它的子樹直到到達葉子節(jié)點,并最終得出結論。

軟組織形變主要是軟組織被提拉或者擠壓,所以軟組織形變的更新主要是AABB樹的更新。發(fā)生形變時,只有接觸點周圍的一些質點發(fā)生了位移形變。相對于整個組織,變形只發(fā)生在局部區(qū)域,所以模型的拓撲結構變化較小,更新時只對發(fā)生了形變的質點進行更新。由于各個節(jié)點形變程度不同,不能用統(tǒng)一函數(shù)來描述。因此,先更新發(fā)生形變處節(jié)點的AABB包圍盒,再對整個包圍盒樹進行更新。父結點的包圍盒更新可通過兩個更新了的子節(jié)點表示。對于整個物體包圍盒的更新,如果檢測到上一層包圍盒沒有相交,其后的包圍盒就無需再更新[16]。

4 虛擬手術仿真結果

使用PC INTEL E3 3.10GHZ,8GBRAM,NVIDIA NVS 300顯卡計算機,以OPENGL圖形庫為基礎,使用VS2010和C++開發(fā)環(huán)境實現(xiàn)軟組織仿真的局部動態(tài)彈性形變。

由于彈簧—質點模型是模擬皮膚表面的形變現(xiàn)象,只是一種近似仿真,而且模型的參數(shù)也無法選取真實軟組織的生物力學參數(shù),所以通過大量實驗,對質點質量、彈簧的彈性系數(shù)和阻尼系數(shù)等進行對比后選取如下:時間步長0.01,質點質量0.01,重力加速度9.8,彈簧的結構彈簧、剪切彈簧和彎曲彈簧的彈性系數(shù)分別為2,2,1.2;阻尼系數(shù)分別為0.4,0.4,0.24。虛擬體彈簧的彈性系數(shù)為表面彈簧彈性系數(shù)的兩倍。

以平面模型為例進行仿真,如圖4所示。分別選取了225個質點,450個面片(a)和625個質點1250個面片(b)的面模型進行試驗。

圖4 不同質點的面模型

以圖4(a)為模型,對固定形變區(qū)域和動態(tài)形變區(qū)域進行對比。當外力較小時,若形變區(qū)域范圍很大,則計算形變花費的時間多,所以我們采用固定形變區(qū)域,如圖5(a)所示;當外力大時,若形變區(qū)域范圍很小,則形變效果較差,所以我們采用動態(tài)形變區(qū)域,如圖5(b)所示。

圖5 形變效果圖

以100個時間步長為例,虛擬一個較小的外力,選用圖5(b)所示模型。在進行形變計算時,選取一個合適的臨界值,頂點的形變大于該值時,激活下一個領域的頂點。若在某個時間步長內,碰撞點在四個方向的形變量均小于該臨界值,則退出循環(huán)。在相同的外力作用下,不同的臨界值不僅影響形變區(qū)域的大小,同時還影響時間步長。臨界值設置越小形變越精細,相應的形變計算步驟越多,時間越長。本文選取臨界值為0.001,虛擬外力為5N的情況下,分別用歐拉法,中點法,四階龍格—庫塔法進行仿真實驗,計算時間如表2所示,仿真效果如圖6所示。

表2 網(wǎng)格模型三種積分方法計算時間比較

實驗結果表明網(wǎng)格質點越多,實時性越差。從3種不同的積分效果去考慮,可以發(fā)現(xiàn),實驗的真實性越來越好,但實時性越來越差。在網(wǎng)格精度相對較高的情況下,歐拉法和中點法能夠達到很好的仿真效果且計算時間能夠達到實時性的需求,但四階龍格—庫塔法的計算時間明顯提高。綜合比較后,在滿足實時性的同時選擇中點法,從而達到較好的仿真精度和效果。

圖6 實驗結果圖

5 結論

本文主要針對虛擬手術仿真過程中實時性和真實性兩方面的需求,采用較能反映軟組織形變“體特征”的質點—彈簧體模型實現(xiàn)了軟組織形變,并在此基礎上提出了更加精確的拓撲結構。為了增強形變過程中的穩(wěn)定性,設置了結構彈簧、剪切彈簧和彎曲彈簧;基于模型形變的動力學方程采用歐拉方法,中點法和四階龍格-庫塔法作為單個質點的運動算法,并通過對虛擬軟組織器官的提拉和按壓手術實驗。結果表明,本文提出的方法能夠很好地滿足模型的精度和仿真實時性要求,增強了形變的仿真體積感,適合于軟組織形變較大時的實時仿真。

[1] 孫連榮,姜元章,任玲.虛擬現(xiàn)實技術及其在高校教學中的應用[J].石油教育,2003(5):41-44.

[2] 戴雯惠.虛擬現(xiàn)實技術的應用現(xiàn)狀及發(fā)展[J].信息技術,2006,35(6):170-174.

[3]Santanu M,Amit R.Simulation of hip fracture in sideways fall using a 3D finite element model of pelvis-femur-softtissuecomplexwithsimplified representation of whole body[J].Medical Engineering&Physics,2007,29:1167-1178.

[4] Zhong Yongmin,Bjian S,et al.Virtual reality simulation of surgery with haptic feedback based on the boundary element method[J].Computers and Structures,2007,85:331-339.

[5] 陳衛(wèi)東,趙成龍.虛擬手術中軟組織形變建模及力反饋算法研究[J].中國生物醫(yī)學工程學報,2013,32(1):113-118.

[6] 王瑞藝.虛擬手術中力反饋的研究和實現(xiàn)[D].鄭州:華北水利大學,2014.

[7]Li Yandong,Zhu Ling.Modeling and simulation of softtissuedeformationbasedonlocaldynamic model[J].Computer Science,2013,40(10):283-288.

[8] Wang P,Becker A A,Jones I A,et al.Virtual reality simulation of surgery with haptic feedback based on the boundary element method[J].Computers and Structures,2007,85:331-339.

[9]馮元禎.生物力學[M].北京:科學出版社,1983:124-158.

[10]Meredith M,Maddock S.Real-time inverse kinematics:the return of the jacobian[J].Department of ComputerScience,TheUniversityofSheffield. 2004:51-60.

[11] 喬冰.基于質點—彈簧模型的人體軟組織形變技術的研究[D].哈爾濱:哈爾濱工程大學,2008.

[12] 褚蓮娣.基于質點一彈簧模型的布模擬方法[J].嘉興學院學報,2002,14(3):57-60.

[13] 潘振寬,高波.手術方針中基于質點一彈簧模型的人體組織變形仿真[J].青島大學學報,2003,18(3):9-14.

[14] 朱元峰,孟軍,謝光華,等.基于復合層次包圍盒的實時碰撞檢測研究[J].系統(tǒng)仿真學報,2008,20(2):372-377.

[15] 鄒益勝,丁國富,許明恒,等.實時碰撞檢測算法綜述[J].計算機應用研究,2008,25(1):8-12.

[16] 趙成龍.虛擬手術中軟組織建模與力反饋算法研究[D].秦皇島:燕山大學,2013.

An Improved Method on Soft Tissue Deformation in Virtual Surgery

HE Wei,WANG Weiping,SHI Weili,MIAO Yu,HE Fei,YANG Huamin
(School of Computer Science and Technology,Changchun University of Science and Technology,Changchun 130022)

Virtual surgery provides the surgery simulation for users on vision and force.It applies computer virtual reality in modern medicine,uses medical image data and reconstructs virtual human soft tissue models.And the simulation and modeling of soft tissue deformation has become one of the most important research contents in virtual surgery system.An improved mass spring model is proposed in this paper for the soft tissue deformation simulation.Bending spring,shear spring and structure spring are added onto the surface model of quadrilateral mesh topological structural for soft tissues,which will make the surface model having volume feature.This model not only has the advantages of traditional spring mass model such as rapid modeling,simple principle and quick simulation,but also has the ability to control the deformation regions,which improves the accuracy of the simulation model,and meets the requirements of the realness and stability of interactive system.Experimental results show that the method is suitable for the real-time simulation of soft tissues with large deformation.It can effectively enhance the volume of deformation simulation and meet the requirements of real-time interactive systems.

virtual reality technology;mass spring model;soft tissue deformation;collision detection

TP317.4

A

1672-9870(2015)06-0118-05

2015-10-30

何巍(1978-),女,博士研究生,講師,E-mail:hw@cust.edu.cn

楊華民(1963-),男,教授,博士生導師,E-mail:yhm@cust.edu.cn

猜你喜歡
結構手術模型
一半模型
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
重要模型『一線三等角』
手術之后
河北畫報(2020年10期)2020-11-26 07:20:50
重尾非線性自回歸模型自加權M-估計的漸近分布
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
論《日出》的結構
3D打印中的模型分割與打包
顱腦損傷手術治療圍手術處理
創(chuàng)新治理結構促進中小企業(yè)持續(xù)成長
主站蜘蛛池模板: 狠狠亚洲五月天| 国产午夜福利在线小视频| 亚洲伊人电影| 久久综合亚洲色一区二区三区| 国产中文在线亚洲精品官网| 四虎国产成人免费观看| 免费在线国产一区二区三区精品| 国产成人亚洲毛片| 九九热这里只有国产精品| 国产97公开成人免费视频| 四虎国产永久在线观看| 免费久久一级欧美特大黄| 99re经典视频在线| 久久婷婷五月综合色一区二区| 一本大道香蕉久中文在线播放 | 日本欧美中文字幕精品亚洲| 国产在线精彩视频二区| 一级毛片不卡片免费观看| 亚洲二区视频| a毛片在线| 国产三级毛片| 麻豆国产在线观看一区二区| 亚洲h视频在线| 日韩欧美在线观看| 国产精品成人久久| 欧美不卡二区| 亚欧乱色视频网站大全| 99精品国产高清一区二区| 国产成人亚洲毛片| 精品中文字幕一区在线| 18禁黄无遮挡免费动漫网站| www.youjizz.com久久| 日韩在线观看网站| 日韩色图区| 被公侵犯人妻少妇一区二区三区| 精品人妻无码区在线视频| 红杏AV在线无码| 国产精品主播| 亚洲系列中文字幕一区二区| 国外欧美一区另类中文字幕| 日本一区二区三区精品视频| 欧美一级视频免费| 日本手机在线视频| 亚洲中文久久精品无玛| 亚洲大学生视频在线播放| 一级毛片高清| 国产精品成人免费综合| 伊人久久青草青青综合| 在线欧美一区| 欧美人在线一区二区三区| 一边摸一边做爽的视频17国产 | 福利视频99| 日韩中文字幕免费在线观看| 亚洲女人在线| 亚洲精品老司机| 在线精品亚洲一区二区古装| 国产免费久久精品99re丫丫一| 久草中文网| 欧美色视频在线| 五月婷婷激情四射| 国产性生大片免费观看性欧美| 中文字幕不卡免费高清视频| 91精品人妻一区二区| 四虎影视8848永久精品| 国产人成乱码视频免费观看| 亚洲无码37.| 日韩精品专区免费无码aⅴ| 欧美成人手机在线视频| 亚洲日本在线免费观看| 美女毛片在线| 免费黄色国产视频| 欧美一区二区三区香蕉视| 99视频国产精品| 久久永久免费人妻精品| 91精品国产自产91精品资源| 天天综合天天综合| 成人小视频在线观看免费| 五月天久久婷婷| 日韩精品成人在线| 亚洲第一色网站| 日韩精品成人网页视频在线| 亚洲国产日韩在线观看|